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I.  INTRODUCTION 


A.  OBJECTIVE 

The  classical  Levinson  algorithm  is  known  to  provide  solutions  to  real  valued,  linear 
systems  involving  Toeplitz  structures.  The  computational  cost  for  these  solutions,  is 
known  to  be  0(k'-).  where  k  indicates  the  filter  order.  It  has  recently  been  proposed  that 
the  classical  algorithm  may  be  transformed  into  2  simpler  algorithms,  using  the  theory 
of  symmetric  polynomials,  and  that  either  of  these  algorithms  can  be  used  to  to  solve  for 
the  predictor  polynomial  of  order  k  at  a  reduced  computational  cost.  (Ref.  1:  p.  47Q| 

These  new  algorithms  are  termed  the  split- Levinson  algorithms  because  their  basis 
is  formed  from  the  concept  of  symmetric  polynomials.  These  are  not  new  in  theory,  but 
the  application  of  the  process  to  linear  prediction  is  a  new  concept.  Symmetric 
polynomials  are  based  on  the  Barlett  Bisection  Theorem  [Ref.  2  •  pp.  1074-1076],  where 
a  system  that  possesses  symmetry  about  a  point,  such  as  a  Toeplitz  matrix,  can  be  de¬ 
composed  into  a  symmetric  and  an  antisymmetric  part.  The  unique  point  of  the  theory 
is  that  either  part  may  be  used  to  solve  the  problem,  or  a  combination  of  both  pans  can 
also  be  useu  in  the  solution.  During  our  rcsca.ch  we  shall  only  consider  real  data  se¬ 
quences. 

The  split-Levinson  case  also  has  a  lattice  structure  as  the  classical  case.  However, 
as  will  be  shown,  the  structure  of  this  lattice  shows  little  resemblance  to  its  classical 
counterpart.  A  derivation  of  a  revised  split  lattice  structure,  and  its  recursive  a'go-ithm 
was  attempted  in  order  to  represent  the  split  lattice  in  a  form  similar  to  that  of  the 
classical  structure.  Although  unsuccessful,  the  derivation  procedure  is  presented  for 
subject  matter  continuity. 

Computer  programs  have  been  written  to  implement  the  new  algorithms  and  com¬ 
pare  them  to  the  classical  algorithms.  Additionally,  computer  programs  are  included  to 
apply  the  split-Levinson  algorithm  for  two  cases,  where  the  computational  efficiency  of 
the  new  algorithm  could  be  of  substantial  benefit.  These  cases  include  the  Moving  Av¬ 
erage  (MA)  problem,  where  the  parameters  of  a  MA  model  must  be  determined  from  the 
given  data,  and  an  extension  of  the  Prony  method  of  spectral  estimation,  where  a  least 
squares  estimation  of  the  presence  of  sinusoids  in  white  noise  is  made  from  the  output 
data  sequence. 
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This  thesis  compares  the  classical  and  lattice  structures  of  the  Levinson  recursion 
formula  given  in  [Ref.  3:  pp.  1 4 5- 1 6T j.  and  examines  not  only  the  formulation  of  the 
recursion  formulas  for  these  algorithms,  but  also  the  complexity  of  the  computations 
and  the  resulting  structure  of  each  of  the  algorithms. 

B.  THESIS  ORGANIZATION 

The  structure  of  the  thesis  is  divided  into  4  chapters,  including  the  Introduction.  In 
Chapter  11  we  will  review  the  classical  Levinson  algorithms.  In  'he  first  case,  the  algo¬ 
rithm  is  obtained  using  the  autocorrelation  function  of  the  input  sequence,  and  in  the 
second  case,  it  is  obtained  using  the  forward  and  backward  error  vectors  of  the  input 
sequence.  In  each  case  we  shall  establish  the  number  of  real  multiplications  required  to 
complete  a  k-th  order  recursion  of  the  respective  algorithm.  As  stated,  the  ultimate  goal 
is  to  establish  the  computational  efficiency  of  the  split-Levinson  algorithm  over  the 
classical  Levinson  algorithm.  Chapter  III  deals  with  the  derivation  of  the  split-Levinson 
algorithms  preceded  by  an  introductory  section  on  symmetric  polynomials.  As  in  Chap¬ 
ter  II.  both  the  autocorrelation  function  and  the  lattice  algorithms  will  be  developed. 
In  addition,  a  comparison  between  the  computational  cost  of  the  Levinson  and  split- 
Levinson  algorithms  and  an  attempt  to  define  the  split  lattice  structure  in  terms  similar 
to  the  Levinson  based  lattice  are  presented. 

In  the  final  chapter  two  practical  applications  of  the  split-Levinson  algorithm  are 
investigated.  These  are:  (1)  the  MA  parameter  estimation  problem,  and  (2)  the  extended 
Prony  method.  In  case  (1).  the  Levinson  recursion  used  to  determine  a  predictor  coeffi¬ 
cient  vector  is  replaced  by  the  split-Levinson  algorithm.  A  comparison  between  the  test 
coefficients  and  the  computed  coefficients  is  presented.  In  case  (2).  an  estimation  of 
sinusoids  in  white  noise  is  performed.  Additionally,  overall  conclusions  of  the  research 
as  well  as  proposed  topics  for  continued  thesis  research  are  presented. 


11.  THE  CLASSICAL  LEVINSON  ALGORITHMS 


The  importance  of  the  Levinson  algorithm  in  linear  prediction  theory  is  well  known. 
The  reason  to  present  the  algorithm  in  us  two  Terms  is  twofold:  th  to  present  certain 
definitions  that  will  be  required  later  in  the  development  of  the  split- Levinson  algo¬ 
rithms.  and  (2)  to  detail  the  computational  complexity  of  the  Levinson  algorithm  lor 
comparison  purposes  to  the  split-Levmson  versions  of  the  algorithm.  In  the  context  of 
our  discussion  within  this  thesis,  we  shall  confine  ourselves  to  the  study  of  autoregressive 


modeling  problems  of  real  sequences  as  in  Figure  1.  [Ref.  3:  p.  152] 


Figure  I.  Autoregressive  Model. 
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We  know  from  linear  prediction  theory  the  augmented  normal  equations  given  by. 

Rx  a  =  rtv  (2.1 1 

can  optimally  be  solved  by  the  Levinson  algorithm,  and  that  this  algorithm  can  be  im¬ 
plemented  with  either  the  autocorrelation  function  or  the  forward  and  backward  pre¬ 
diction  error  vectors  of  the  input  sequence  [Ref.  3:  pp.  152-1  ~0j . 

A.  THE  LEVINSON  ALGORITHM 

In  order  to  examine  the  computational  complexity  of  the  Levinson  recursion,  it  is 
necessary  to  formulate  the  recursive  algorithm,  and  to  determine  the  number  of  real 
multiplications  and  additions  required  to  complete  the  algorithm.  First,  construct  a 
Toeplitz  matrix  from  the  sequence  s(t).  of  length  N.  defined  as  R.  =  0  <  /,  j  <  A], 

where  the  elements  of  the  matrix  are  the  autocorrelation  lags  given  by  [Ref.  2:  p.  646] 


R, 


v- ;  -1 


.=.) 


(2.2) 


then,  the  predictor  vector  a  can  be  determined  as  a  solution  to  the  system  defined  by  the 
matrix  equations 

C/?J[aJ  =  [^.0,0 . ,0]r  (2.3) 

where  a,  is  the  prediction  error  norm,  and  is  defined  as 


ak  = 


k 

*  +  £  aki  R-l 

i=l 


(2.4) 


It  is  recalled  from  linear  prediction  theory,  that  given  a  positive-definite  matrix  R. 
of  order  k+  1,  the  kth  order  a  coefficient  vector  can  be  computed  recursively  from  the 
nested  Toeplitz  submatrtces,  and  their  respective  successive  predictor  vectors,  a.  The 
well-known  Levinson  recursion  formula  is  this  solution,  and  has  the  form 

^kl  =  ^k—\,i  Pkl*k—\  .k—i  /— 0. 1  A'  (-.5) 


with  the  conditions  that  at0  =  1.  and  at_,*  =  0.  The  parameters,  p*  =  ut.»,  are  called  re¬ 
flection  coefficients,  also  PARCOR  coefficients,  because  they  represent  the  partial  cor- 
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relation  between  the  zero-th  and  the  k-th  clement  of  the  prediction  vector  with  the  effect 
of  all  the  intermediate  elements  removed. [  Ref.  4:  p.  53| 

To  construct  the  Levinson  recursion  we  must  use  the  prediction  error  norm  re¬ 
lationship 

ak  =  1 1  ~ 

and  the  identity 

A*— 1 

a  ;  pk  =  -  2^  Rk-flk- u 

i=0 

to  define  the  recursion  variables.  Consider  the  following  definition  as  it  applies  to  the 
Levinson  recursion  [Ref.  1:  p.  472j. 


*-i 

/ *  =  ~  X;  Rk-iak-\,i 

;=0 

=  Ojc-lPk 

and  solving  for  pk  from  Eq.<2.8),  we  have 


(2.8) 


(2.4) 


The  error  norm  ak  can  be  written  in  terms  of  the  the  normalizing  term  /,  by  rewriting 
L'q.  (2.6),  and  making  a  substitution  from  Eq.  (2.9) 

<**  =  (!-  p\)°k-\ 

-pk{ak_xpk)  12.10) 

=  ak- 1  —  Pk*-k 

Combining  Eqs.  (2.5),  (2.8),  (2.9),  and  (2.10),  we  have  the  basis  for  the  Levinson  algo¬ 
rithm.  and  it  is  summarized  in  Table  2  of  Appendix  A. 

B.  LEVINSON  LATTICE  REALIZATION 

If  we  are  given  a  real  sequence  of  signal  values  s(0)  s(  1) . s(N-l).  and  it  is  known 

that  s(t)  =  0.  for  -!  >  t  and  i  >  .V.  then  for  the  linear  prediction  problem  of  order  k  we 
find  it  necessary  to  find  a  set  of  real  numbers  a^.  ... ,  ak,  that  will  minimize  the  for¬ 

ward  and  backward  prediction  errorvectors  using  a  linear  combination  of  the  past  signal 
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vectors,  [f  we  call  the  forward  prediction  error  vector fX‘)  and  the  backward  error  vector 
bjr).  and  define  them  in  terms  of  the  a.,  coefficients.  [Ref.  2:  p.  646 J  we  have 

jr 

Ur)  =  V a  sir  -  i)  1 2.1 1 ) 


bkir)  —  2  akj<-i  ~  i)  (2.12l 

then  it  turns  out  that  the  same  real  numbers.  ak,  .  will  provide  the  solution  to  either  of 
the  forward  or  backward  prediction  problems,  ( i.e. .  minimize  the  squared  Euclidean 
norm  of  both  J]  and  b,). 

Let  a.  be  defined  as  the  squared  norm,  that  is 

From  Eqs.  (2.11)  and  (2.12)  forming  the  first  three  terms  of  each  error  vector  we  have 
the  following. 


M 0)  =  rf*oilO)  +  ak]s(  -1)  +  ...  +  akks(  -k) 

M I)  =  «*o5(l)  +  a*,s(0)  +  ak:s(  - 1 )  +  ....  +  akks(  1  -  k) 
fk(2)  =  ak0s{ 2)  +  a*,s(l)  +  ak2s{ 0)  +  ...  +  akks( 2  -  k) 


(2.15) 


t>k{0)  =  akks(°)  +  akJ<- 15(  -*)  +  <Zk*-A  “2)  +  -  +  %35(  -*) 
**(1)  =  +  ak.k-\s(° )  +  ak.k- 2S(  _1)  +■•■  +  ajwill  -  *) 

bk(2)  =  a**s(2)  +  akJc_,s(  1)  +  akJc-2s(0)  +  ...  +  akQs{ 2  -  k) 


(2.  IS) 
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If  we  examine  the  elements  of  these  two  vectors  we  can  see  they  are  related  in  that 
each  can  be  derived  from  the  other  by  reversing  the  order  of  the  a  coefficients.  If  we  form 
the  Euclidean  norm  of  each  vector,  ;!f„!  and  b.  .  we  see  that  the  4.-th  predictor  vector 
[a*,-,  ih. . it»Jr  minimizes  the  error  norm,  and  =  \l\". 

From  (Ref.  3:  pp.  156-157],  we  can  use  the  Levinson  algorithm  to  deline  the 
recursion  formula  for  the  forward  and  backward  prediction  errors  given  by 

Jl(r)  +  pA-i('  -  1  >•  ,  , 

If  we  let  the  following  definition  apply  to  the  lattice  version  of  the  Levinson  algo¬ 
rithm 


v+jr-2 

'■k  =  °k-\Pk~-  Y,  /*_i(O&*_i(0  (2-21) 

f=l 

then  using  Eqs.  (2.9),  (2.10),  (2.20),  and  (2.21)  we  can  summarize  the  Levinson  lattice 
algorithm  in  Table  3  of  Appendix  A. 

Even  though  the  lattice  algorithm  is  implemented  directly  from  the  data  samples,  its 
computer  implementation  will  be  more  complex  because  of  the  vector  manipulations 
that  must  occur  in  each  iteration.  The  lattice  structure  defined  by  Eq.  (2.20)  is  shown  in 
Figure  2. 

In  summary,  we  discussed  the  Levinson  algorithm  which  uses  the  autocorrelation 
elements  in  its  recursion  and  the  related  lattice  structure  which  uses  the  input  data  di¬ 
rectly  in  its  formulation.  In  terms  of  the  computational  complexity,  both  algorithms  re¬ 
quire  real  multiplications  of  the  order  k2.  as  detailed  in  Tables  2  and  3  of  Appendix  A, 
in  order  to  realize  a  k’*  order  predictor  filter. 
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III.  THE  SPLIT-LEV fNSON  ALGORITHMS 


The  split-Levinson  algorithms  are  based  on  the  theory  of  symmetric  and  antisym¬ 
metric  polynomials.  We  know  that  for  an  k-th  order  real  autoregressive  process,  the 
normal  equations  are 


where  B  represents  one-half  of  the  vector  components  of  a*,  and  B'  represents  the  re¬ 
versal  of  the  vector  components  of  B.  Using  Eqs.  (3.3)  and  (3.4).  we  can  transform  the 
normal  equations  into 
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Raf  =  C^/2,0.0. ...  a;2]r 
Rif  =  [<x/2.U.U. . 

Therefore,  we  can  see  some  favorable  consequences  of  these  revised  normal 
equations,  and  their  solutions.  First,  either  the  symmetric  or  antisymmetric  form  will 
give  the  same  solution,  and  second,  because  of  the  symmetry  of  the  predictor  vectors, 
we  need  only  solve  for  one-half  of  the  predictor  coefficients. 

Similar  to  the  Levinson  algorithm  we  now  proceed  to  develop  the  split- Levinson 
algorithms  from  the  input  sequence  autocorrelation  function  and  the  predictor  error 
vectors. 

A.  SPLIT-LEVINSON  ALGORITHM 

The  predictor  polynomial  a4(z)  is  defined  as 

k 

ak(:)  =  YjakiZ~‘  (3-6) 

i=0 

relative  to  the  given  Toeplitz  matrix  of  autocorrelation  lags.  Denote  the  reverse  of  our 
predictor  polynomial  as  a[r)  =  r-*a*("'),  and  the  predictor  polynomial  has  been  shown 
to  obey  the  recursion  [Ref.  3:  pp.  156-157] 

a*(;)  =  «*_i(2)  +  /V_'a(2)  (3.7) 

and  the  reverse  polynomial  of  Eq.  (3.7)  is 

ak{:)  =  z~]ak_l{z)  + pj.a^iz)  (3.8) 

We  now  want  to  form  a  new  polynomial  from  the  given  predictor  polynomial  that 
will  form  the  basis  of  the  split- Levinson  algorithm.  It  is  desired  to  show  that  the  deter¬ 
mination  of  the  coefficients  of  this  polynomial  will  allow  us  to  recover  the  original  pre¬ 
dictor  polynomial,  and  at  the  same  time  be  more  computationally  efficient.  We  define 
the  symmetric  polynomial  as  Pt(z),  and  the  antisymmetric  polynomial  as  P[c)(z)  ,  and  we 
desire  them  to  be  of  the  form  [Ref.  1:  p.  472] 
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fy-)  =  Pk;= 


P*-' 


Recall  from  Eq.  (3.4)  and  (3.5)  that  the  symmetric  and  antisymmetric  predictor  cocfTi- 
cients  are  composed  of  two  vectors  that  are  reverses  of  each  other,  and  we  will  define 
these  vectors  so  that  they  obey  the  relationships 


Pki  Pkjc—i 

l.7>  id) 

Pk.i  Pk.k—i 


(3.10) 


Consider  the  mathematical  interpretation  of  making  the  autocorrelation  matrix.  Rt 
a  singular  matrix.  If  the  reflection  coefficient  pk  is  made  ±  1.  then  this  corresponds  to 
an  element  of  R,  making  the  matrix  singular.  For  this  reason  we  shall  designate  the 
symmetric  and  antisymmetric  predictor  polynomials  as  singular  predictor  polynomials 
[Ref.  2:  p.  472]  and  from  Eq.  (3.9)  they  are  defined  as 

/*(*') +  1(->  -  ... 

M  >  (-'-ll) 

Pk  U)  =  ^-iU)  -  Z  rtfc_i(-> 


Also,  these  singular  predictor  polynomials  are  seif-reciprocal  [Ref.  2:  p.  472]  because  of 
their  symmetry  and  may  be  expressed  in  the  following  forms 

w-'V'i  l3P) 

From  Eq.  (3.1 1)  we  have 

2~'4-i(2)  =  W~<W-) 

If  we  add  Eqs.  (3.7)  and  (3.8).  and  make  a  substitution  from  Eq.  (3.13)  we  have 

<*£*)  +  «*(*)  =  *“'«*-!(*)  +  Pk“k-  |U)  +  ak-\(z)  +  Piflk-M) 

=  (1  +  pk)ak^{z)  +  (1  +  pk)z~'ak_\{z)  (314) 

=  f kP *(-) 
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where  we  have  defined  /„  as 


'■k  ~  1  +  Pk  (3.15) 

In  a  similar  fashion  we  can  solve  for  the  antisymmetric  normalized  singular  predictor 
poly  normal  by  subtracting  Eqs.  (3.7)  and  (3.8).  and  substituting  from  Eq.  (3.13)  we  have 

«*(-')  -  u,(e)  =  ^P^(z)  (3.16) 


where 


•<<*)_  i 

*  Pk 


(3.17) 


Similar  to  the  predictor  polynomial  a»(z),  we  can  define  the  singular  predictor  coef¬ 
ficient  vectors  for  Eq.  (3.9)  as  [Ref.  2  :  p.  472] 


Pa  “  LPkO<  Pki . Pkk\T 

Ja)  W-\T 
Pa  —LPkO-Pkl . PkkJ 


(3.18) 


Since  we  want  the  split- Levinson  normal  equations  to  be  of  the  form 


Pk*k  =  l>*.0.0, ...  ,0]r 

(3.19) 

R*a*  =  [0,0,...,(7jr 

(3.20) 

then  from  Eq.  (3.14)  and  (3.16)  the  singular  predictor  polynomials  are 


pk(z)  = 


■+■  ■ 


P(k\Z)  = 


gfcU) 

;(a) 


Ak 

to. 

:(<J) 


(3.21) 


Since  ak(z)  is  a  polynomial  formed  from  the  predictor  coefficient  vector  that  is  a  solution 
to  Eq.  (2.3),  it  follows  that  Pk{z)  and  P\°\z)  are  solutions  to  the  Toeplitz  system  described 
by  Eqs.  (3.19)  and  (3.20).  [Ref.  1:  p.  472] 
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KiPk  =  R' 


r,?:'  =  r 


a(. 

f  k 


a;. 

(j) 


(3.221 


Normalizing  Eqs.  (3.19)  and  (3.20)  by  /,  and  /.'f,  the  split- Levinson  normal  equations  in 
matrix  form  are  expressed  as 

^i:Pk=  [>*.0.0,  ...  ,Tjr 

. U‘'' 

where  we  have  defined  the  modified  prediction  error  norm  [Ref.  2:  p.  651] 


(3.24) 


If  we  expand  the  matrix  expressions  in  Eq.  (3.23).  the  modified  error  norms  may  be  ex¬ 
pressed  as  finite  sums  of  the  predictor  coefficients 


Tk 


T 


(a) 

k 


(3.25) 


where  R,  is  the  i-th  autocorrelation  element  of  the  k-th  sub  matrix. 

Since  the  symmetric  and  antisymmetric  polynomials  are  closely  related,  we  shall 
derive  only  the  symmetric  polynomial  recursion  equations,  and  then  simply  present  the 
results  for  the  antisymmetric  case. 

The  final  step  in  the  derivation  is  to  derive  a  three  term  recursion  formula  for  the 
symmetric  polynomial.  From  Eq.  (3.11)  and  (3.14)  we  have  the  surprising  result  that  the 
predictor  polynomial  a„(r)  can  be  obtained  from  a  linear  combination  of  successive  sin¬ 
gular  predictor  polynomials  [Ref  2.  p.  472].  First,  form  /V,(r)  from  Eq.  (3.1 1),  and  then 
eliminate  a„(:)  using  Eq.  (3.14) 
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(3.2m 


PM{z)  =  flfeU)  +  -  ak(z) 

=  ak(: I  +  r-'[/*^(-)  ~  a*U.i] 

=  <  1  -  z_1)a*U)  +  :~'/.kPki:) 

(1  -  .-■')«*(:)  =  ^+1(--)  -  z~‘kkPk(z) 

If  we  replace  k  by  k-1  in  Eq.  (3.26)  above  we  also  have 

We  now  form  our  recursion  formula  by  mutiplying  Eq.  (3.1 1)  by  (1  —  z1 ),  and  use  Eqs. 
i  3.13).  (3.26).  and  (3.2")  to  eliminate  p,  and  all  at  predictor  polynomials. 

(1  -  :~‘)ak{:)  =  (1  -  z_i)ak_,(z)  +  p*( I  -  z~')z~' ak_x(z) 

=  (1— Z  *  ;(.**)  "b  Pk{  ^  Z  ~ 

=  (1  -  z'’')«k_,(r)  +  pk[Pk(z)  -  ak_x(z)  -  z~]  Pk(z)  4-  z_1ak_,(z)]  (3.2S) 
=  «*_i(-)[l  -  z"‘  -Pk  +  Z~'pk ] 

=  «*_!<-)[(  I  ~  Z-‘)(l  -  Pk)] 

If  we  now  substitute  for  (1  — z'1)u»(z'1),  (1  —  z-'K-ilz’1)  from  Eqs.  (3.27)  and  (3.28).  and 
eliminate  p»  using  Eq.  (3.15).  we  can  complete  the  derivation 

Pk+lC)  -  z-]AkPkiz)  =  lPk(z)  -  z",/k_,^_1(z)](2  -  /.*)  +  [ z.kPk(z)  -  Pk(z)](l  -  z"1) 

=  2 Pk(z)  -  Ak(z)Pk(z)  -  2z~]/.k_xP *_,(z)  +  z"  Vk-A-iU) 

+  /^(z)  -  AW  -  z-'s.kPk(z)  +  z->,(z)  (3.29) 

^+1(z)  =  (1  -  z~l)Pk(z)  +  z~,/*_lP*_1(z)[/k  -  2] 

=  (!  +  z~])Pk(z)  -  «kz~' Pk_x(z) 

Taking  the  inverse  Z  transform  of  Eq.  (3.29),  we  have  the  three  term  recursion  formula 
for  the  singular  predictor  coefficients 

Pki  =  Pkt  +  PkM-i  ~  *kPk-\jc-i  (-3  :,°> 

where  the  recursion  parameter  a*  is  defined  as 

«*-/*-.[2->-*]  (3-31^ 

We  note  that  r*  is  determined  from  Pk(z)  from  Eq.  (3.23),  and  therefore  we  conclude  that 
the  singular  predictor  polynomials  can  be  recursively  computed  from  Eq.  (3.30).  How- 
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ever,  the  recursion  parameter  a*  is  not  quite  in  the  correct  form.  From  Eqs.  i2.1<0, 
(3.15),  and  (3.31)  we  can  alternatively  compute  a,  from  [Ref.  2:  p.  473] 


= 


13.32) 


The  dual  relationships  for  the  antisymmetric  split-Levinson  formulas  can  be  derived 
by  following  a  procedure  similar  to  the  one  presented  above.  It  suffices  to  replace  the 
quantities  /„  r.a.. by  their  antisymmetric  duals,  i.e.,  pvf,  and  use  the  following  anti¬ 
symmetric  initial  conditions.  (Ref.  2:  p.  649] 


A,o  =  0 

<j)  i 

P  10  =  1 

(J)  , 

p  ii  =  -  i 


Tu-> 

”0 


=  & 


(3.33) 


Recursive  equations  for  the  symmetric  split-Levinson  algorithm  are  summarized  in 
Table  4  of  Appendix  A.  Examining  the  entries  in  Tablp  4,  we  see  that  a  full  iteration 
loop  of  the  algorithm  requires  approximately  t  real  multiplications.  However,  because 
of  the  symmetry  of  the  singular  predictor  coefficients,  we  only  have  to  perform  one-half 
of  these  calculations.  Therefore,  for  a  k-th  order  filter  we  need  to  make  on  the  order  of 
Ay  2  real  multiplications.  The  <)  function  in  Table  4  is  used  to  distinguish  between  even 
and  odd  orders  of  the  indexing  variable. 

The  FORTRAN  program  SPLIT,  in  Appendix  B,  estimates  the  predictor  coefficients 
using  the  Levinson  and  split-Levinson  algorithms.  Figure  3  is  a  graphical  comparison 
between  the  known  test  filter  coefficients  of  SPLIT,  shown  by  the  solid  curve,  and  the 
filter  coefficients  computed  by  the  Levinson  and  split-Levinson  algorithms,  shown  by  the 
dashed  curve.  We  now  undertake  the  derivation  of  the  lattice  form  of  the  split-Levinson 
formulas  to  verify  the  numerical  complexity  of  that  method,  and  to  investigate  the 
symmetric  and  antisymmetric  lattice  structures  compared  to  the  Levinson  lattice  forms. 

B.  SPLIT  LATTICE  ALGORITHM 

We  begin  the  split  lattice  derivation  by  introducing  the  symmetric  and  antisymmetric 
error  vectors  x*.  xf  [Ref.  2:  p.  648].  If  we  use  previously  established  singularity  concepts, 
and  substitute  ±  l  for  p„  in  Eq.  (2.20)  for  the  symmetric  and  antisymmetric  error  vectors, 
x*  and  Xf\  respectively,  we  have 
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o 

CJ 


Figure  3.  Levinson  /  Split  Levinson  Coefficient  Comparison. 

=/*-,(') |('-D  (3  3J) 

4flV)  =A_,(o  -  Vi('~  » 

As  in  the  split- Levinson  case,  we  shall  proceed  with  the  derivation  of  the  symmetric 
split  lattice,  and  present  the  antisymmetric  lattice  results  at  the  end  of  the  derivation 
with  any  significant  changes  noted. 
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If  we  extend  the  singular  polynomial  concept  to  the  singular  predictor  coefficients, 
we  can  start  with  the  Levinson  coefficient  recursion  formula  and  substitute  -  1  for  p. 

11  k:  =  ak- 1  i  +  Pkak-  .k-i  ,3-3 5 1 


and  substituting  for  pv 

Pki  =  +  «*-!.*-/  <1361 

If  we  write  Eq.  (3.34)  representing  the  time  index  (t)  with  the  subscript  (i)  we  have  an 
algorithm  that  is  more  easily  adapted  for  computers. 

%=A-u  +  i,H^  l3-37' 

Now,  comparing  Eqs.  (3.36)  and  (3.37)  we  have  a  direct  correlation  between  the  two, 
and  from  the  split-Levinson  equation  for  the  forward  error  vector,  we  can  write  the  dual 
split  lattice  equation  for  .*•*(/) 


k 

xk{t)  =  X  pkl  sir  -  i)  (3.38) 

;=0 

Since  Eq.  (3.3S)  is  in  the  form  of  a  convolution  sum  we  can  apply  Z  transform  theory 
to  see  if  any  inferences  can  be  made 


A3(c) 
Xk(z) 
5( -*) 


Pk(:)S(z) 

Pki-) 


(3.39) 


From  Eq.  (3.39)  we  can  conclude  that  the  symmetric  polynomial  constitutes  the  Z 
transform  of  the  transfer  function  formed  from  the  error  vector  and  input  sequence  z 
polynomials.  Now  we  can  use  the  previously  derived  split-Levinson  algorithm,  and  in¬ 
verse  transform  it  to  obtain  the  lattice  error  vector  recursion  algorithm. 

Repeating  the  split-Levinson  recursiun  we  have 

Pk+l{:)  -  (1  +  2~')/V=)  -  *kz~'pk- i(z)  ,  4()) 

=  Pki2)  +  2~]pki2)  ~  ^k2'‘ Pk-\iz) 

Multiplying  Eq.  (3.40)  by  Slz)  and  using  Eq.  (3.39)  for  P*(z)S(z)  we  have 
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(3.41) 


S(;)Pk+l 


A’, 


fc+l 


(-)  =  5(r)[Pk(r)  4-  r  1  Pk{:)  -  a^k-i'-)] 

u)  =  A’*u)  +  r]xk{z)  -  xkr'1  xk_.u) 


Applying  the  inverse  Z  transform  to  each  side  of  Eq.  (3.41)  we  have  the  singular  pre¬ 
dictor  error  vector  recursion  formula  (Ref.  2:  p.  650 J 

.v^ilr)  =  .vt.(/l  +  r-l.r..(/)  -  xkz~\xk_.{t)  (3.42) 


The  symmetric  lattice  structure  given  by  Eq.  (3.42)  is  shown  by  Figure  4.  (Ref.  2: 
p.  65b] 

Erom  Eq.  (3.32)  we  know  that  the  recursion  parameter  a*  is  defined  as  and 

since  i.  appears  in  the  recursion  formula  for  the  singular  error  vector,  we  need  to  solve 
for  it.  We  begin  with  the  Levinson  error  norm  equation 


-°k- 


**■(*  ~  p\)°k- i 

°k  =  -<**-,(  1  +  P*)(l  -  Pk) 

=  -  Pk) 


1  +  pk 

( 1  -  Pk)  *  2*k 


(3.43) 


where  we  have  substituted  r*  from  Eq.  (3.32),  and  2cr*_;(  1  -  pk)  is  defined  as  ||jr*;i2  [ Ref. 
2:  p.  650]. 

The  initial  conditions  for  Eq.  (3.42)  must  be  examined  because  most  cases  are  trivial 
except  for  the  case  of  k  =  0.  From  Eq.  (3.38)  we  have 

*o(0-Poo*M.  (3-44) 


and  from  [Ref.  2:  p.  648]  we  define  px  -  2. 

The  antisymmetric  duals  are  very  similar  to  the  symmetric  case,  and  can  be  formed 
by  replacing  the  symmetric  variable  by  its  antisymmetric  counterpart,  i.e.,  xp(t)  for 
xt(i),  -  p„  for  p»,  etc.  The  initial  conditions  for  the  antisymmetric  case  are  given  below 


4%)  =  o 

X(k%)  =  s(i)  -  s{t  -  1) 


(3.45) 


The  antisymmetric  lattice  structure  given  by  the  antisymmetric  dual  of  Eq.  (3.42)  is 
shown  in  Figure  5.  The  split-Levinson  lattice  formulas  given  by  Eqs.  (3.37),  (3.43),  and 
(3.44)  are  summarized  in  Table  5  of  Appendix  A.  If  we  examine  the  number  of  real 
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to 


Figur:  4.  Symmetric  Lattice  Structure. 
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multiplications  given  in  Table  4  and  Table  5,  then  we  can  deduce  the  following  conclu¬ 
sions.  The  split  versions  of  the  Levinson  algorithm  do  produce  a  reduction  in  the  com¬ 
plexity  of  the  calculations  when  compared  to  the  classical  versions.  The  split-Levinson 
produces  a  reduction  of  one-half,  and  the  Levinson  produces  a  reduction  of  3  2  [Ref.  2: 
p.  645].  The  FORTRAN  program  SLATIS.  Appendix  C,  implements  the  Levinson  and 
split-Levinson  lattice  algorithms,  and  a  graphical  comparison  between  the  known  test 
coefficients,  shown  by  the  solid  curve,  the  coefficients  estimated  by  the  Levinson  lattice 
algorithm,  shown  by  the  dashed  curve,  and  the  coefficients  estimated  by  the  symmetric 
split-Levinson  algorithm,  shown  by  the  dotted  curve,  is  presented  in  Figure  6. 

The  split  lattice  structures  shown  in  Figure  4  on  page  19,  and  in  Figure  5  on  page 
20  show  that  the  classical  lattice  structure  appears  to  be  lost  in  the  new  algorithm.  The 
distinct  advantage  of  the  original  lattice  structure  is  the  modularity  of  the  filter.  In  order 
to  retrieve  this  appealing  feature  of  the  lattice  filter  we  shall  now  proceed  to  derive  a 
revised  version  of  the  split  lattice  structure,  and  see  if  it  can  have  a  form  similar  to  the 
classical  structure. 

C.  SPLIT  LATTICE  REVISED  STRUCTURE 

To  begin,  consider  the  second  order  classical  lattice  structure  derived  from 
Figure  2  on  page  8.  We  can  write  the  following  equations  for  the  first  and  second  stage 
forward  and  backward  prediction  errors, 

/;(«)  =  s(n)  +  pls(n  -  1) 

£,(/?)  =  p,S(/»)  +  5(/l  -  n 

/;(«)  =/i(«)  +  PigM  ~  1)  (  '  ’ 

=  p/i(«)  +  Si(«-  1) 


Now  substituting  the  equations  for/fn)  and  forg,(n)  into  the  equations  for/2(n)  and 
£T:(n).  solving  for  the  second  stage  forward  and  backward  prediction  filter  errors  and 
taking  the  Z  transforms  yields  the  transfer  functions 


F^z) 

~S{z) 


=  1+2  '(p,  +  p,p2)  +  P2Z 


(3.47) 


and. 


G2(z) 

5(2) 


=  Pi  +Z  ’(Pi  +  P\Pi)  +  2 


(3.48) 
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Figure  6.  Levinson  vs.  split-Levinson  Coefficient  Comparison, 
which  obviously  yield  the  second  order  predictor  transfer  function  A2(z)  and  its  reverse 
version  A2(z),  where  a20  =  1,  aa  =  p2.  «  Pi  -  P\Pi  -  PiO  -  Pi)-  'Sow  fon™11?  t*16  ^ird 

order  symmetric  polynomial  P2(z)  from  our  second  order  example,  we  have  [Ref.  1:  p. 


If  we  define  a0  =  1.  and  i<;,  4-  a: )  =  p,  .  then 


P-J:)  =  1  t  p.:  -  P\Z  '  - 

=  1  ^  /'.J-'  -r-  --“'I;-. 


.-3 

-1 , 


(3.5<U 


Define  the  following  terms: 


Fi(r)  =  1  +  pt: 

G,(r)  =  ;_l(l  */»,;) 


(3.51) 


therefore, 

F3(e)  =  F,(e)  +  r-2C,(r)  (3.52) 


Equation  (3.52)  defines  the  revised  symmetric  split-lattice  structure,  and  Figure  7  gives 
a  graphical  representation  of  that  structure. 


Fi(2) 

_ ^ _ 1 

x(k)^ 

-  1 
7 

-  2 

'w/cyry— i 

z 

Figure  7.  Proposed  Split  Lattice  Configuration 


In  order  to  show  that  the  preceding  classical  structure  can  be  equated  with  the 
symmetric  polynomial  derived  earlier,  it  is  now  necessary  to  form  the  forward  and 
backward  prediction  error  transfer  function  of  the  new  split  lattice  structure  and  com¬ 
pare  them  to  the  Z  transform  of  the  symmetric  polynomial.  Therefore,  from  Figure  7 
and  Eq.  (3.52),  we  can  write  P}(z)  as 

P3(z)  =  1  +  A’,*"1  +  z"2(A',  +  z~ ')  (3.53) 
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Now,  if  we  compare  Eq.  (3.50)  to  £q.(3.53),  we  see  that  A",  =  p,.  But,  from  Eq.  (3.40) 
we  know  that  p{  =  (a,  +-  a2),  and  substituting  for  a,  and  a2  from  Eq.  (3.4S)  and  (3.49) 
yields 


*'i  =  P\  =Pi  +P2<1  -Pi)  (3.54) 

Let  us  now  redenve  the  symmetric  polynomial  from  the  revised  split  lattice  structure. 
From  Eq.  (3.50)  and  (3.54)  we  have 

P3(j)  =  1  +  (p,  -  p2  -  +  z~2[(p,  -  Pi  -  pxp2)  4-  (3.55) 

Since  we  have  found  that  the  symmetric  lattice  can  be  restructured  to  a  form  similar 
to  the  classical  lattice,  the  next  logical  step  is  to  find  a  recurrence  relation  for  the  new 
lattice.  Let  us  consider  a  5*  order  symmetric  polynomial  to  determine  the  step-down 
recursion  procedure,  given  by 

/^(r)  =  1  +  p5l2~!  +  pS2z~2  +  pS2z~2  +  p5lz^  +  z~s  (3.56) 

From  Eq.  (3.51)  we  have 

f2(z)  =  1  +  p5,c  '  +  p$iz  2  |3  5~) 

G2(z)  =  z~2  +p512-1  +p$2 

As  per  the  observations  made  in  Figure  7  on  page  23  and  Eq.  (3.51),  we  have  the  lattice 
reflection  coefficient  for  the  second  stage,  K2  =  pS2 ■  Now  reduce  the  order  of  FLi)  to  find 
the  first  stage  reflection  coefficient  using  the  standard  inverse  Levinson  recursion  [Ref. 
4:  pp.  156-157], 


fiU) « 


F2(z)  -  K2G2(z) 


■  1  + 


i  -  Kt 

Pri 

1  +  A2 


-i 


therefore 


A,= 


P  51 

1  +  A*2 


(3.58) 


(3.59) 


Now  rewriting  the  equations  for  F{{z)  and  F2(z)  using  the  derived  reflection  coefficients, 
we  have 
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(3.601 


F](z)  =  l  +  A[Z 

<J|U)  =  -  '  +  A, 

F2{z)  =1  +  Aj(  1  -t-  A2)z  +  KiZ 
G2{z  )  =  z  +  Aj(l-t-A2)z  +  A; 


and  we  know  that 

/>5U)  =  f2(z>  +  z*3G;(r) 

=  1  +  A,(l  +  A;)z_1  +  A\z'2  +  A2z-3  (3.61) 

+  A";(  l  4-  A"2)z  +  z 

Equating  the  terms  in  Eq.  (3.56)  and  (3.61),  we  have  p,A  =  A,(l  +  A',)  and  p!2  =  A\. 
Knowing  the  values  of  A',  and  A:,  we  now  can  form  a  two  stage  symmetric  lattice  similar 
to  that  in  Figure  7  on  page  23  to  implement  P,(z)  .  However,  we  are  interested  to  find 
if  we  can  recursively  obtain  the  lower  orders  P4(z).  Ps(z).  etc.  or  the  higher  orders 
P,(z).  P-{z).  etc.  from  Ps  (z).  In  an  attempt  to  form  PJz),  we  use  the  standard  forward 
recursion  [Ref.  3:  pp.  156-15”)  to  obtain 

F2(z)  «  f2iz)  +  K\z~'G2(z) 

-  1  +  A’,(  1  +  A2)z-1  +  A\z-2  +  A2A3z-1  (3.62) 

+  A  VA3(  l  +  A2)z  +  Aj  z 


G3(z)  =  i  3  +  A, (l  t-  AV)z 

+  K2z~'  +  A2A3;':  (3.63) 

+  A,  A3(  1  +  A2)z  +  A3 

Fornung  the  symmetric  polynomial  P„(z)  from  Eqs.  (3.62)  and  (3.63)  we  have 

P6(2)  =  F3(2)  +  z'3G3(z) 

=  1  +(A|  +  A]  At  +  A2A3)z 

+  (K2  +  A,  Kj  +  A,  A2A3)z-2  +  A3z-3  (3.64) 

+  A3z  3  +  (A2  +  A,  A3  4-  A[ A2A3)z 
+  (A,  +  A,A2  +  A2A3)z-5  +  z~6 


Comparing  Eq.  (3.64)  to  the  symmetric  form  of  P6(z),  we  have 


Pei  =  (^i  +  ^1  i-  A.\ k\ ) 


P?\  PsiPn 

1  +  Pi2  If  7^52 


1 

+  T  PszPei 


(3.65) 


The  one-half  enters  into  Eq.  (3.65)  because  we  know  that  for  even  orders  the  polynomial 
coefficients  are  symmetric  about  the  center  element,  and  they  must  be  shared  in  the 
matrix  equations.  We  shall  now  expand  Eq.  (3.65)  to  attempt  to  develop  a  recursive 
algorithm  for  the  sixth  order  coefficients  from  the  fifth  order  coefficients.  Expanding  Eq. 
(3.65)  we  have 


(1  +  Pn)Pti  =Ps  1  +P$\P52  + 


1  +P52 
2 


jl>52  +  Ptl  -  *sPtll 


(3.66) 


where  the  term  in  brackets  is  an  expansion  of  the  coefficient  recursion  formula,  Eq. 
(3.30),  for  p6i.  Collecting  terms  we  have 

(1  +Psi)P6\  =/>5l(l  +/>52)-(l  +  Pil)  T  P^lPsi  +  Psl  ~ 

*5  (3-67) 

Pei  ■/»st  +  /»jjOs:  - 
Substituting  for  p!2  from  Eq.  (3.30) 

P*1~P*1  +  Pa\  -«uP3i  (3.6S) 

From  Eq.  (3.68)  we  can  observe  that  the  a*  recursion  parameter  and  the  number 
of  previous  coefficients  required  to  be  known  are  increasing  in  order,  and  it  appears  that 
a  simple  recursive  algorithm  based  on  the  above  approach  is  not  possible.  Note  that 
although  the  new  lattice  structure  does  not  appear  to  be  order  recursive,  we  can  express 
a  given  order  symmetric  or  antisymmetric  lattice  structure  in  a  more  conventional  form. 

In  summary,  we  know  that  the  split-Levinson  algorithm  is  a  viable  replacement  for 
the  classical  algorithm  because  of  its  computational  efficiency.  We  have  studied  both 
autocorrelation  and  data  (or  lattice)  realizations  of  the  split-Levinson  algorithm.  An  at¬ 
tempt  to  derive  a  recursive  split  lattice  algorithm  yielded  a  classical-like  lattice  structure, 
but  it  is  not  recursive  in  order.  Further  investigation  is  necessary  in  this  direction.  We 
now  need  to  test  the  split-Levinson  algorithms  on  some  signal  processing  applications. 
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IV.  APPLICATIONS  OF  THE  SPLIT-LEVINSON  ALGORITHM 


In  this  chapter,  we  apply  the  split-Levinson  algorithm  in  (1)  the  MA  parameter  es¬ 
timation  problem,  and  (2)  the  extended  Prony  method  of  spectral  line  estimation.  Before 
we  take  up  these  two  applications  we  examine  the  algorithm's  usefulness  if  the  original 
filter  has  coefficient  sy:.  metry.  i.e. ,  the  impulse  response  of  a  linear  phase  FIR  filter. 

A.  HANKEL  AND  TOEPLITZ  MATRICES 

In  previous  derivations  we  ’nave  assumed  the  FIR  filter  equation  to  be  non- 
symmetric.  Let  us  now  investigate  the  problem  where  the  filter  equation  is  symmetric, 
i.e.,  of  the  form 


y(n)  =  s{n)  +  a{s(n  -  1)  +  a2s(n  -  2)  + 
+  1)  +  aks(n  —  k) 


(4.1) 


By  definition,  a  symmetric  polynomial  is  self-reciprocal,  that  is 

tf/cU)  =  (J-) 

Therefore,  from  the  Levinson  algorithm,  predictor  polynomials  are  known  to  obey  the 
recurrence  relation 


ak(:)  =  <?*_,(;)  4-  pkz  'ak_{{z)  (4.3) 

and  in  our  special  case  we  have 

ak(z)  =  a*_,(r)  +  pkz~]ak_{(z)  =  a{z)  J( 

=  (1  +  pkz~])ak_^z) 

In  order  to  formulate  a  set  of  equations  similar  to  the  split-Levinson.  it  is  necessary 
to  derive  the  normal  equations  for  our  special  case,  and  compare  them  to  the  standard 
equations,  in  order  to  develop  the  recursive  algorithms.  Since  the  predictor  coefficients 
in  a  recursive  algorithm  produce  estimates  of  s(n)  as  the  algorithm  steps  through  its  re¬ 
cursive  steps,  denote  this  estimate  as  i(«).  In  vector  form,  we  then  have 
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5(rt|rt  —  1)  =  —  [s(rt  -  1 )  ^(/T  —  2 )  ...  5(0)] 


a\ 

a2 

«3 

«2 

a. 


(4.5) 


To  derive  the  normal  equations  we  must  find  the  minimum  mean  squared  error  (MSE) 
from  the  equation  for  the  error, 

e{n)  =  s(r.)  —  s{n[n  —  1)  (4. 6' 

The  minimum  mean  squared  error  is  found  by  squaring  the  error  term,  and  then  differ¬ 
entiating  the  squared  term  with  respect  to  the  given  ak  vector.  Combining  these  two  ev¬ 
olutions  we  have  the  following  equations, 

MSE  =  J 

=  £!>>)]  (4.7) 

=  £[(5(«)  -  (5(«)|/I  -  1))J] 

To  obtain  the  normal  equations,  we  carry  out  the  following  steps: 

77  -  0  =  -2£[s(«)s*_,]  -2£Is[_,sa_,] 

C  a 

E[s{k) sk_,s[_,]a  -  £[5(n)s*_,]  |J-8» 

rf-') 

From  the  split-Levinson  recursion  formulas  we  know  that  the  singular  symmetric 
polynomial,  /^(r),  is  defmed  as  the  following. 

Pk(:)  =  a^iz)  +  z~]  a^iz)  (4.9) 

Since  our  predictor  polynomial  is  symmetric,  it  is  a  reasonable  question  to  ask  if  sym¬ 
metric  polynomials  also  obey  this  recursive  relationship.  It  is  noted  as  an  immediate 
consequence  of  the  recursive  problem,  all  of  the  preceding  polynomials  will  also  be 
symmetric.  Therefore,  we  check  the  singular  predictor  polynomial  recursion  to  see  if  it 
holds  when  the  original  polynomial  is  itself  symmetric 
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a*(‘)  =  <**-!(*)  +2  4_lW 

=  a*_i(2)  +  z_1z_\_,(r) 

=  1  +  c^z-1  +  a2z-x  +  ■•■  +  z~{k~])  +  ;~k[l  +  tJ,z_1  +  u2z_~  4-  ... 
=  l+  (l+CJ])z  +(<2|+£!j)Z  ~  4-  ...  +  z 


-k  HI") 
:  A] 


Now  that  we  see  that  the  Levinson  recursive  equation  holds  for  a  symmetric 
polynomial,  we  derive  the  recursive  relationships  for  our  polynomial  using  what  is 
known  from  the  split-Levinson  equations.  We  have  defined  the  symmetric  polynomial 
P.(z)  to  be  a  normalized  combination  of  a  non  symmetric  polynomial,  at{ z),  and  its  re¬ 
ciprocal  image.  at(z)  in  the  form  of, 

'kPkiz)  =  pH’)  +  4u)  HI  l) 

By  direct  substitution  it  is  a  trivial  matter  to  show  that  this  relationship  also  holds  for 
a  symmetric  polynomial,  at(z). 

In  order  to  develop  the  recursion  for  the  symmetric  polynomial,  it  is  necessary  to 
express  the  desired  linear  predictor  in  terms  of  the  previous  two  predictors.  To  this  end 
use  Eq.  (4.10),  to  form  a,.,(r),  and  substitute  from  Eq.  (4.11)  to  perform  this  task. 

+  2-'ak-i(z) 

=  <**(2)  +  2-1£*(r) 

=  ak(z)  +  z~'[/kak(z)  -  ad.-)] 

=  ak{:){  1  -  z_1)  +  :~'/.kak{z) 

Solving  for  (1  -  z*')at(z),  and  forming  the  quantity  (1  -  z-')a*_,(z)  we  have 

(1  —  z~l)ak(z)  =  ak+1(z)  —  z~l  /.kak(z) 

-1  -1  HI  3) 

(1-2  )«k-l(2)  =  flfc(z)  ~  Z 

These  relationships  will  now  allow  us  to  form  the  three  term  recursion  for  the  given 
symmetric  polynomial  from  Eqs.  (4.3),  (4.1 1),  (4.12),  and  (4.13) 

«k(?)  =  ak- 1(2)  +  pA-  1(2)  ,  .  . 

<Wr)  =  (1  +  z  )ak(z)-<xkz  afc_,(z) 


where  we  have  defined  a* 


aJe  —  *k— 1(2  4) 
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From  Eq.  (4.14)  we  can  see  that  the  coefficient  recursion  formula  is  the  same  form 
as  Eq.  (3.29),  and  we  can  deduce  that  the  split-Levinson  algorithms  will  work  equally 
well  for  symmetric  polynomials  that  describe  unknown  filters  as  it  does  for  polynomials 
that  are  not  symmetric. 

B.  FIR  MOVING  AVERAGE  PARAMETER  ESTIMATION 

If  we  consider  an  FIR  filter  with  an  input  sequence  given  by 
sr=  [.s(«)s(«  -  l)...s(«  -  m)],  and  an  output  y(n)  given  by 

1 1 

v(rt)  =  Vu,  s{n  -  i)  (4.16) 

<=o 

then  we  can  develop  the  necessary  equations  to  estimate  the  moving  average  parameters, 
and  solve  for  the  FIR  filter  coefficients.  The  algorithm  to  estimate  the  predictor  coeffi¬ 
cients  can  be  defined  as  follows: 

Let  the  three  predictions,  y?(n),  s?{n),  and  sp(rt),  represent  the  m-th  order  predictions 
of  the  forward  estimate  of  y,  the  forward  estimate  of  s,  and  the  backward  estimate  of  s 
respectively.  [Ref.  6:  p.  1  ] 


yf'{n)  =  'Yj>is{n- i)  (4.1") 

(=0 

m 

sf  (n)  =  s(n  —  i)  1 4. 1 S) 

1=0 

m 

s™(n  —  m)  =:'S^dl  s(n  —  m  +  i),  (4.19) 

1=0 

where  the  coefficient  vectors  are  defined  as. 

br  =  [1  —  bQ  -bx 

cf=  [0  1  -c,  ... -cj  (4.20) 

dr=  [0-  cm  -  cm_,  ...  -  c,  ] 

Without  going  through  all  the  details  for  the  MA  parmeter  recursions,  we  can 
summarize  the  recursion  formulas  for  the  predictor  coefficients  as  [Ref.  6:  p.  2] 
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Notice  that  the  predictor  vector  d"  is  not  included  in  the  preceding  definitions  since  it  is 
the  reverse  c" 

If  we  examine  Eq.  (4.21)  we  see  that  the  recursive  relationship  for  c"  is  a  statement 
of  the  Levinson  recursion,  since  X",  and  Av"  are  the  m-th  order  reflection  coefficients. 
Therefore,  we  can  apply  the  split- Levinson  algorithm  to  solve  for  c".  form  d"  .  and 
recursively  determine  b”.  finally,  from  the  theory  of  Moving  Average  processes,  bm  =  a’ 
The  LORTRAN  program  MAV1.  in  Appendix  D,  uses  a  25-th  order  LIR  test  filter 
to  to  generate  a  test  sequence,  and  the  results  are  given  in  Table  6  of  Appendix  A. 

figure  S  is  a  graphical  comparison  between  the  known  test  coefficients,  shown  by 
the  solid  curve,  and  the  computed  filter  coefficients,  shown  by  the  dashed  curve.  The 
curves  are  fitted  to  the  linear  magnitudes  of  the  coefficients  by  interpolating  spline 
techniques. 

C.  EXTENDED  PRONY  METHOD 

The  estimation  of  the  existence  of  sinusoids  in  the  presence  of  noise  is  a  common 
occurrence  in  signal  processing  applications.  In  this  simulation,  we  will  show  that  the 
spilt- Levinson  algorithm  can  be  directly  implemented  in  the  process  to  determine  the 
approximate  frequencies.  The  noise  is  zero  mean,  unit  variance,  and  white  in  nature. 

In  this  application  of  the  split-Levinson  algorithm,  we  attempt  to  approximately  fit 
p  exponentials  to  a  data  set  of  N  samples.  We  have  the  constraint  that  .V  >  2 p.  and  a 
least  squares  estimation  procedure  is  used.  We  begin  by  defining  the  estimate  of  our  set 
of  data  samples.  [Ref.  7:  p,  1404] 


p 

x  =  Ybmznm  n  =  0.....Y  —  1  (4.22) 

m=l 

where  bm  =  Am  exp(jdJ/2,  and  zm  =  exp(j2nfmAt).  The  zm  s  are  roots  of  uni*  modulus  with 
arbitrary  frequency  and  occur  in  complex  conjugate  pairs  as  long  as  /„  ¥=  0  or  l/2Ar 
Therefore,  in  order  to  solve  for  the  p  sinusoids,  we  must  solve  for  the  roots  of  the  fol¬ 
lowing  polynomial.  [Ref.  7:  p.  1406] 


31 


Figure  8.  MA  Coefficient  Comparison. 

2r 

T(r)  =  V  akzlp~k  =  0  (4.23) 

*=0 

The  roots  can  be  of  unit  modulus,  and  occur  in  complex  conjugate  conjugate  pairs  if 
we  constrain  the  polynomial  coefficients  to  be  symmetric  about  the  center  element  of  the 
polynomial  [Ref.  7:  p.  1407],  This  is  an  exact  ocurrence  in  the  symmetric  and  anti  sym¬ 
metric  polynomials  of  the  split-Levinson  algorithm,  as  long  as  the  order  of  the 
polynomial  is  even,  that  is 


ir  = 


1  o. 


up-\ 


(4. 24) 
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Note  that  the  last  element  of  the  coefficient  vector  is  a.  2  rather  than  a.  because  of  the 
symmetry  of  the  polynomial,  and  that  symmetric  polynomials  only  guarantee  that  if  a 
root  r  occurs,  then  so  does  its  reciprocal  e;1 2  [Ref.  7:  p.  140’’]. 

The  program  EPRONY1,  Appendix  E.  utilizes  the  split-Levinson  algorithm  to  ap¬ 
proximate  the  p  order  sinusoids  to  the  given  data  set.  A  summary  of  the  simulation  cases 
studied  are  presented  in  Table  1,  and  the  graphical  results  follow. 


Table  1.  SUMMARY  OF  TEST  CASES 


Case 

Number 

Constants 

Variables 

i 

Npts 

Test  frequencies 

fs 

SNR(-10,  0,  10  dB) 

-> 

Test  frequencies 

f 

J  s 

SNR 

NPTS 

3 

Test  frequencies 

fs 

SNR 

Npts 

Filter  Order 

4 

Test  frequencies 

fs 

Npts 

SNR 

Filter  Order 

SNR  (-10,0.10  dB) 

1.  Simulation  Parameter  Definitions 

All  simulation  cases  are  done  in  the  presence  of  white  noise,  and  a  minimum 
possible  separation  frequency  for  the  input  sinusoids  was  determined.  All  plots  are 
sinusoid  magnitude,  in  a  linear  scale,  versus  digital  radian  frequency,  6  for  0  <  6  <  n. 

2.  Simulation  Results 

We  begin  the  spectral  line  estimation  simulations  with  all  parameters  fixed,  and 
then  selectively  choose  a  parameter  to  vary  and  observe  the  effects  of  these  changing 
parameters.  We  begin  with  two  sinusoids  of  f  —  50  Hz,  f  —  75  Hz,  which  yield 
t>,  =  1.396  radians,  and  =  2.0944  radians,  respectively.  The  number  of  data  points 
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(NPTS)  is  set  at  1500,  and  the  filter  order  (M)  is  chosen  to  be  4  indicating  the  presence 
of  two  sinusoids.  Now,  we  chose  to  vary  the  signal-to-noise  ratio  (SNR)  for  10  dB.  0  dB, 
and  -10  dB,  shown  in  Figure  9. 

Figures  9  (a)  and  9  (b)  show  that  lowering  the  SNR  from  10  dB,  Figure  9  (a), 
to  0  dB.  Figure  9  (b),  causes  the  estimation  to  worsen,  and  both  indicate  low  frequency 
estimation  error.  From  both  of  these  cases  we  can  deduce  the  presence  of  two  sinusoids, 
but  in  Figure  9  (c)  it  appears  that  only  one  sinusoid  is  present,  and  the  low  SNR  has 
caused  spectral  estimation  to  fail.  The  conclusion  for  the  first  case  is  that,  the  better  the 
SNR,  the  better  the  spectral  estimate. 

For  the  second  case  we  selected  NPTS  as  the  variable  parameter,  held  the  filter 
order  and  sinusoid  frequencies  constant  as  before,  and  set  the  SNR  at  0  dB.  From  Fig¬ 
ures  10(a),  (b),  and  (c)  we  clearly  see  that  the  better  estimation  occurs  with  NPTS  = 
1000.  Figure  10(b),  because  of  the  equal  amplitude  and  accurate  frequency  estimation 
as  compared  to  Figures  10  (a)  and  10  (c).  All  plots  show  low  frequency  error,  but  also 
suggests  that  simulations  should  be  conducted  with  NPTS  set  to  1000-1500  points  for 
the  best  results.  This  case  provides  the  rationale  for  the  value  of  NPTS  for  all  other 
simulations. 

In  the  third  simulation  all  parameters  are  fixed  as  in  the  second  case,  and  M  is 
varied  for  4,8.  and  10.  From  Figures  11(a),  (b),  and  (c)  we  can  see  that  although  there 
are  only  two  sinusoids  present,  the  estimation  plots  show  M  2  sinusoids  for  all  values 
of  filter  order.  Each  plot  has  frequency  components  in  the  vicinity  of  the  actual  fre¬ 
quencies,  but  they  also  give  false  indications  of  spectral  lines.  If  we  were  making  a  deci¬ 
sion  on  the  number  of  frequencies  based  on  the  magnitude  plots  of  Figures  11(b)  and 
(o,  then  signifigant  errors  would  be  introduced.  In  this  case  it  is  obvious  that  unless  the 
exact  number  of  sinusoids  present  is  known,  then  the  estimation  technique  will  fail. 

In  the  final  simulation  we  introduce  two  additional  sinusoids,  and  examine  the 
effects  of  SNR  on  spectral  line  estimation.  The  constants  for  this  case  are  f  —  35  Hz. 
[  =  85  Hz,/}=  105  Hz,/4  =  175  Hz,  NPTS  =  1500,  M  =  4,  and  525  Hz.  Digital 
frequencies  are  0.419,  1.010,  1.496,  and  2.094  radians  per  second  respectively.  As  in  case 
l,  SNR  takes  on  the  values  of  10  dB,  0  dB,  and  -10  dB.  In  Figure  12(a),  where  the  SNR 
is  10  dB,  we  get  good  results  with  near  uniform  amplitude  estimation,  and  digital  fre¬ 
quencies  that  are  close  for  all  frequencies.  However,  from  Figures  12  (b)  and  12  (c),  we 
can  see  that  the  spectral  line  for /,  is  missing,  and  an  errant  line  appears  at  approximately 
2.8  radians.  Both  of  the  figures  show  unequal  amplitude  indications,  but  3  of  the  4 
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spectral  lines  are  in  close  proximity  to  their  actual  values.  Troni  this  case  we  draw  the 
conclusion  that  the  spectral  line  estimation  performance  deteriorates  at  low  S\Rs. 


Figure  9.  Spectral  Estimation:  Filter  Order  =  4;  Data  Record  Length  =  1500; 
SNRs:  (a)  10  dB,  (b)  0  dB.  (c)  -10  dB. 
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Figure  10.  Spectral  Estimation:  Filter  Order  -  4:  SNR  -  0  dB,  Data  Record 
Lengths:  (a)  500),  (b)  1000,  (c)  3000), 
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Since  the  overall  purpose  of  the  simulation  was  to  test  the  applicability  of  the 
split- Levinson  algorithm  to  the  test  cases,  then  it  has  been  shown  that  the  split-Levinson 
will  produce  estimates  for  the  respective  cases.  However,  if  we  examine  the  accuracy  of 
the  low  signal-to-noise  ratio  cases,  we  see  that  the  new  algorithm  suffers  a  similar  fate 
as  the  classical  case  in  that  it  is  difficult  to  accurately  estimate  the  correct  sinusoidal 
frequencies  in  the  presence  of  noise. 

D.  CONCLUSIONS 

The  split- Levinson  algorithm  has  been  shown  to  be  computationally  more  efficient 
than  its  classical  counterpart.  We  can  say  that  the  application  of  the  split-Levinson  al¬ 
gorithms  to  practical  applications  in  lieu  of  the  Levinson  algorithm  can  be  advanta¬ 
geous,  and  the  computational  cost  can  be  reduced  significantly  for  large  order  systems. 
Additionally,  the  split-Levinson  algorithms  are  applicable  to  problems  where  we  want 
to  model  MA  parameters,  and  perform  spectral  estimation  using  the  Prony  method. 

We  note  the  following  restrictive  areas  for  the  new  algorithm  that  could  make  it 
unsuitable  for  certain  signal  processing  applications. 

1.  Non-recursive  split-lattice  algorithm. 

2.  Computational  accuracy  degradation  when  performing  spectral  estimation  in  low 

signal-to-noise  cases. 

3.  Complexity  of  symmetric  and  antisymmetric  lattice  structures. 

Since  we  have  shown  that  the  split-Levinson  algorithm  is  a  viable  substitute  for  the 
Levinson  recursion,  it  is  reasonable  to  consider  areas  of  this  topic  for  further  research. 
We  know  that  the  symmetric  lattice  structure  can  be  expressed  in  a  classical  form  for 
given  filter  orders,  however,  a  recursive  algorithm  for  this  new  structure  has  been  elusive. 
We  propose  that  the  existing  recursive  algorithm  for  the  singular  predictor  coefficients 
should  be  studied  to  see  if  a  redefinition  of  equation  parameters  can  extract  a  new  re¬ 
cursive  algorithm  for  the  new  symmetric  lattice.  Additionally,  the  algorithm's  poor  per¬ 
formance  in  low  signal-to-noise  ratio  test  cases  of  the  extended  Prony  method  is  similar 
to  the  performance  of  the  classical  algorithm.  Therefore,  techniques  to  improve  the 
classical  algorithm's  performance,  as  detailed  in  [Refs.  8.9],  may  be  investigated  for  ad¬ 
aptation  to  the  split-Levinson  algorithm. 
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APPENnix  A.  TABULAR  SUMMARY  OF  ALGORITHMS 


The  tables  given  in  this  Appendix  were  taken  f-om  [Ref.  2:  pp. 648-674),  and  are 
presented  for  the  convenience  of  the  reader. 


[  Table  2.  THE  LEVINSON  ALGORITHM 
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Table  3.  THE  LEVINSON  LATTICE  ALCORITHM 
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Table  4.  THE  SPLIT-LEVINSON  ALGORITHM 
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Table  5.  THE  SPLIT  LATTICE  ALGORITHM 
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Table  6.  MOVING  AVERAGE  TEST  RESULTS 
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APPENDIX  B.  SPLIT-LE VINSON  PROGRAMS 

PROGRAM  TO  CALCULATE  THE  NTH  ORDER  FIR  PREDICTOR  FILTER  USING 
THE  SYMMETRIC  AND  ASYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIALS, 
THE  SPLIT- LEVINSON  RECURSION  FORMULA,  AND  THE  AUTOCORRELATION 
FUNCTION  OF  THE  INPUT  SEQUENCE. 

VARIABLE  DEFINITIONS 

SIGMAN  -  N-TH  DEGREE  NORM  OF  THE  FILTER. 

KEVEN  -  INTEGER  VARIABLE  USED  TO  CONTROL  ACCESS  TO 

SUBROUTINE  EVEN  WHEN  THE  INDEX  VARIABLE  K  IS  AN 
EVEN  INTEGER. 

LAMDA  -  REAL  VECTOR  USED  WHEN  DEFINING  THE  SINGULAR 

LAMDAS  PREDICTOR  POLYNOMIAL  (PK(Z))  IN  TERMS  OF  THE 

NORMALIZED  SYMMETRIC  AND  ANTISYMMETRIC  PARTS  OF 
THE  PREDICTOR  POLYNOMIAL  AK(Z). 

LAMDA(K)  =  1  +  RHO(K) 

LAMDAS(K)  =  1.  -  RHOS(K) 

ALPHA  -  REAL  VECTOR  USED  TO  SIMPLIFY  THE  THREE-TERM 

ALPHAS  RECURRENCE  RELATION  FOR  THE  SINGULAR  PREDICTOR 
POLYNOMIALS  OF  THE  SAME  TYPE. 

ALPHA(K)  =  LAMDA(K-1)*( 2.  -  LAMDA(K) ) ,  OR 
ALPHAOO  =  TAU(K)/TAU(K-1) 

ALPHAS(K)  =  TAUS(K)/TAUS(K-1) 

KODD  -  INTEGER  VARIABLE  USED  TO  CONTROL  ACCESS  TO  THE 
SUBROUTINE  ODD  WHEN  THE  INDEX  VARIABLc  K  IS  AN 
ODD  INTEGER. 

TAU  -  REAL  VECTOR  OF  "MODIFIED  NORM  VALUES".  THE 

TAUS  VALUES  ARE  CALCULATED  FROM  A  SUMMATION  OF 

PRODUCT  TERMS  OF  THE  AUTOCORRELATION  LAGS,  AND 
THE  COEFFICIENTS  OF  THE  SINGULAR  PREDICTOR 
POLYNOMIALS. 

RHO  -  REAL  VECTOR  OF  REFLECTION  COEFFICIENTS  RHO(l), 

RHOS  RHO(2) ,. .  .  ,RHO(N). 

N  -  DESIRED  ORDER  OF  THE  PREDICTOR  POLYNOMIAL. 

C  -  REAL  VECTOR  OF  AUTOCORRELATION  LAGS  R(0),R(1), 

R(2) .  ,R(N) 

P  -  ARRAY  OF  SINGULAR  PREDICTOR  POLYNOMIAL 

PS  COEFFICIENTS  FROM  P(0,0) . . P(N+1,N). 

A  -  ARRAY  OF  PREDICTOR  POLYNOMIAL  COEFFICIENTS. 

AS  EACH  I-TH  ROW  OF  THE  ARRAY  CONTAINS  THE 

PREDICTOR  POLYNOMIAL  COEFFICIENTS  FOR  THE  I-TH 
ORDER  PREDICTOR  POLYNOMIAL. 

T  -  INTEGER  VARIABLE  USED  IN  THE  SUBROUTINE  ODD 

IN  THE  COMPUTATION  OF  THE  TAU(K)'S. 

L  -  DUMMY  VARIABLE  USED  DURING  THE  CALCULATION 

OF  THE  SYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIAL 
COEFFICIENTS. 

VARIABLE  DECLARATIONS 
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DEAL  C( 0: 100) ,P(0: 100,0: 100) ,TAU(0: 100) 

REAL  A( 0: 100,0: 100) , LAMDA! 0: 100) ,RH0! 100) .SIGMAN 
REAL  PS( 0: 100,0: 100),AS(0: 100,0: 100) 

REAL  RHOS( 100), ALPHAS! 100 ),TAUS!0:  100), SIGMAS 
REAL  AR( 0:  30) ,W(0:  5000) ,S(0: 5000) ,SIGMA(0:  100) 

REAL  R( 0:  100) , ALPHA! 100) , LAMDAS! 0:  100) 

REAL  AA! 0:  100,0:  100) , GAM( 50), LAMK , LAM! 100) 

INTEGER  M , LL , IX , T , KODD , KEVEN , L , N 
OPEN’! UNIT=4,BLANK='  ZERO'  ) 

C  INITIALIZE  FILTER  ORDER 
READ! A, 100 )N 
LL  =  N 
IX  =  1913 
M  =  3000 
WRITE! 6 , 300) 

WRITE! 6 , 400  jN 
WRITE! 6 , 450  )M 
DO  1  I =0 , M 

CALL  GAUSS! IX, 1.  ,0.  ,V) 

W(  I )  =  v 
1  CONTINUE 

CALL  INPUT! LL, W , AR , S ,M) 

C  INITIALIZE  AUTOCORRELATION  VECTOR  C 
CALL  ACORR! C , S ,N+1 ,M) 

WRITE! 6, 2 100) 

C  INITIAL  CONDITIONS  FOR  THE  SYMMETRIC  AND  ASYMMETRIC  PREDICTOR 
C  POLYNOMIAL  CALCULATIONS. 

P(0,0)  =  2. 

P(l.l)  *  I- 
TAU(O)  =  C(0) 

LAMDAIO)  =  1. 

KODD  =  1 
KEVEN  =  2 
A(N,0)  =  1.0 
PS!0,0)  =  o. 

PS! 1,1)  =  -1. 

TAUS(O)  =  C!0) 

LAMDAS! 0)  =  1. 

AS(N,0)  =  1.  0 

CALL  LEV!C, GAM, N.AA, LAM, SIGMA) 

WRITE! 6, 1700) 

WRITE (6, 1800) 

DO  20  J=1,N 

WRITE! 6, 1900)N,J, SIGMA! J) ,LAM!J) ,AA(N,J) 

20  CONTINUE 

C  SYMMETRIC  &  ASYMMETRIC  SPLIT-LEVINSON  RECURSION 
WRITE! 6 , 800) 

WRITE! 6, 850) 

WRITE! 6, 900) 

DO  99  K=1,N 
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p(K,cn  =  l. 

PS ( K ,  0)  =  1. 

C  CALL  STATEMENTS  FOR  EVEN  OR  ODD  VALUES  OF  INDEX  K 
IF( K.  EQ.  KODD)THEN 

CALL  AODD( C , PS , K , N.TAUS ,T) 

CALL  ODD(C,P,K,N,TAU,KODD,T) 

ELSEIF( K. EQ. KEVEN)THEN 
CALL  AEVEN( C , PS , K ,N ,TAUS ,T) 

CALL  EVEN( C , P , K , N , TAU , KE VEN , T) 

END  IF 

ALPHA( K)  =  TAU(K)/TAU(K-1) 

ALPHAS! K)  =  TAUS(K)/TAUS(K-1) 

C  LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  40  1=1, T 

P( K+l , I )  =  P( K , I )  +  P!K,I-1)  -  ALPHA(K)*P(K-1,I-1) 

PSIK+l.I)  =  PS(K,I)  +  PS(K.I-l)  -  ALPHAS! K)*PS(K- 1,1*1) 

C  DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTGR  COEFFICIENTS 
IF(  I.  LT.  T  .  OR.  I.  EQ.  K)GOTO  40 
L  =  K+l 
DO  50  J=I+1,K 
P!L,J)  =  P(L,L-J) 

50  PS( L , J)  =  -PS(L,L-J) 

40  CONTINUE 

LAMDA(K)  =  2.  -  (ALPHA(K)/LAMDA!K-1)) 

LAMDAS(K)  =  2.  -  ( ALPHAS! K)/LAMDAS!K-1)) 

C  REFLECTION  COEFFICIENT  CALCULATION 
RHO(K)  =  LAMDA(K)  -  1. 

RHOS(K)  =  1.  -  LAMDAS!K) 

WRITE! 6, 1000)K,RHO!K) ,RHOS(K) ,GAM(K) 

99  CONTINUE 

C  CALCULATION  OF  N-TH  ORDER  NORM  ISIGMAN)  AND  N-TH  ORDER  PREDICTOR 
C  COEFFICIENTS,  A(N, 1) ,A!N,2) ,. . . ,A(N,N) 

SIGMAN  =  LAMDA! N ) *TAU ( N ) 

SIGMAS  =  LAMD AS (N)*TAUS(N) 

WRITE! 6 , 1 100) 

WRITE! 6, 600) 

WRITE! 6 , 1200) 

DO  60  1=1, N 

A!N,I)  =  ACN.I-l)  +  PCN+1,1)  -  LAMDA!N)*P!N,I-1) 

AS!N,I)  =  AS(N, 1-1)  +  PS!N+1,I)  -  LAMDAS(N)*PS(N , 1*1) 

WRITE! 6 , 1300) I ,TAU! I ) ,TAUS( I ) , ALPHA! I ) , ALPHAS! I) ,P( 1+1 , I ) , 

+PS( 1+1,1) ,A(N,I) , AS( N , I ) 

60  CONTINUE 

100  FORMAT! 12) 

200  F0RMAT!F12. 6) 

300  FORMAT! '1') 

400  FORMAT! ’  FILTER  ORDER  =  ’,13) 

450  FORMAT!'-’,’  NUMBER  OF  SAMPLE  POINTS  =  ’,15) 

600  FORMAT! ,103X,' FILTER  COEFFICIENTS') 

700  FORMAT! 5X, 13, 1 IX, F10.  4,21X,F10.  4) 
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C  J  CM 


800  FORMAT( ,2 IX, 'REFLECTION  COEFFICIENTS') 

850  FORMAT! ,  25X, '  SPLIT-LEVINSON'  ) 

900  FORMAT! 5X, ' INDEX' , 8X ,' SYMMETRIC ' , 9X ,' ANTISYMMETRIC ' ,9X, 

+' LEVINSON' ) 

1000  FORMAT! 5X, 13, 10X.F12.  6,12X,F12. 6.12X.F12. 6) 

1100  FORMATC ,5X, 'FILTER  PARAMETERS  FROM  SPLIT  LEVINSON  RECURSION') 

1200  FORMATC  5X, ' INDEX'  ,6X, 'TAU(K)' ,4X,'TAU*(K) '  ,  4X , ' ALPHA! K) ' ,4X, 

+  ' ALPHA*( K) ' ,4X, ’PCX)' ,4X,'P*(K)' ,6X, ' SYMMETRIC' , 6X , ' ASYMMETRIC '  ) 
1300  FORMATC 5X, 13, 6X.F12.  6,3X,F12.  6,3X,F12.  6,4X,F12.  6 , 3X.F12.  6 , 2X, 
+F12.  6 , 2X ,F10.  4 , 5X,F10. 4) 

1700  FORMAT! '-' ,5X, 'FILTER  PARAMETERS  FROM  LEVINSON  RECURSION') 

1800  FORMAT!'-' ,8x,' INDEX' ,8X, ' SIGMA(K) ' ,5X, ' LAMDA(K) ' ,8X,'FILTER 
+COEFFICIENTS ' ) 

1900  FORMAT! 8X, 12, 13, 7X.F12.  6.6X.F12.  6.12X.F12.  6) 

2000  FORMAT! 5X, ' SPLIT-LEVINSON  RECURSION  CALCULATIONS') 

2100  FORMAT! 5X, ' INDEX’ ,5X,' KNOWN  COEFFICIENTS' , 5X , ' AUTOCORRELATION 
+FUNCTION  C ! K ) ' ) 

200  FORMAT! 2X, 13, 4X.F12.  6) 

300  FORMAT! 5X , 13 ,40X,F12.  6) 

WRITE! 6 ,300) 

END 

SUBROUTINE  AODD! C , PS ,K ,N ,TAUS ,T) 

C  THIS  SUBROUTINE  CALCULATES  THE  ANTISYMMETRIC  "MODIFIED 

C  NORMS"  WHEN  THE  INDEX  K  IS  AN  ODD  INTEGER. 

REAL  C ( 0 : 100) ,PS(0: 100,0: 100) ,TAUS(0: 100) 

INTEGER  T 
T  =  iK-lJ/2 
TEMP  =  0. 

DO  5  1=0, T 

5  TEMP  =  TEMP  +  ! C( I )  -  C(K-I ) )*PS(K, I) 

TAUS( K)  =  TEMP 

RETURN 

END 

SUBROUTINE  ODD! C , P , K ,N ,TAU , KODD ,T) 

C  THIS  SUBROUTINE  CALCULATES  THE  "MODIFIED  NORMS"  (  TAU!K)'S) 

C  WHEN  THE  INDEX  K  IS  AN  ODD  INTEGER. 

REAL  C ! 0 :  100), P! 0:  100,0:  100),TAU!0:  100) 

INTEGER  T 
T  =  (K+l)/2 
TEMP  =  0. 

DO  15  1=0, T-l 

15  TEMP  =  TEMP  +  (C(I)  +  C( K- I ) )*P! K, I ) 

TAU(K)  =  TEMP 
KODD  =  KODD  +  2 
RETURN 
END 

SUBROUTINE  AEVEN( C , PS , K , N.TAUS ,T) 

C  SUBROUTINE  CALCULATES  THE  VALUE  OF  THE  ANTISYMMETRIC 

C  "MODIFIED  NORMS"  (TAUSCK)'S)  WHEN  THE  INDEX  K  IS  AN  EVEN 

C  INTEGER. 
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REAL  C( 0:  100) ,PS(0:  100,0:  100) ,TAUS(0:  100) 

INTEGER  T 
T  =  K/2 
TEMP=  0. 

PS(K,T)  =  0. 

DO  25  1=0, T-l 

25  TEMP  =  TEMP  +  (C(I)  -  C(K-I) )*PS(K, I) 

TAl’S(K)  =  TEMP 

RETURN 

END 

SUBROUTINE  EVEN( C,P,K,N,TAU, KEVEN ,T) 

C  SUBROUTINE  CALCULATES  THE  VALUE  OF  THE  "MODIFIED  NORMS" 

C  (TAU(K)'S)  WHEN  THE  INDEX  K  IS  AN  EVEN  INTEGER. 

REAL  C( 0: 100) ,P(0: 100,0: 100) ,TAU(0: 100) 

INTEGER  T 
T  =  K/2 
TEMP=  0. 

DO  35  1=0, T-l 

35  TEMP  =  TEMP  +  (C(I)  +  C(K-I' )*P(K, I) 

TAU(K)  =  TEMP  +  C(T)*P(K,T; 

KEVEN  =  KEVEN  +  2 

RETURN 

END 

SUBROUTINE  INPUT( LL ,W , AR ,S ,M) 

C  SUBROUTINE  TO  GENERATE  THE  INPUT  SEQUENCE  FROM  A  GIVEN  FIR 

C  FILTER  AND  ZERO  MEAN,  UNIT  VARIANCE  WHITE  NOISE. 

REAL  W( 0: M) ,AR( 0: LL) ,S( 0:  M) 

C  CALCULATE  INPUT  SEQUENCE  VALUES  BETWEEN  0  AND  FILTER  ORDER. 

C  TEMP  =  0. 

S(0)  =  W( 0) 

DO  45  K=1,M 

S(K)  =  W(K)  -0. 6*W( K- 1 )  +  0. 8*S(K-1) 

45  CONTINUE 
RETURN 
END 

SUBROUTINE  ACORR( C , S , NLAG , M ) 

C  SUBROUTINE  TO  CALCULATE  THE  AUTO  CORRELATION  FUNCTION  OF  THE 
C  GIVEN  INPUT  SEQUENCE. 

REAL  C( 0: 100), S(0: 5000) 

INTEGER  T 
DO  5  1=0, NLAG 
TEMP  =  0. 

DO  15  T=0,M-1-I 

TEMP  =  TEMP  +  S(T)*S(T+I) 

15  CONTINUE 

C( I)  =  TEMP*( 1.  /FLOAT  (M-l-I)) 

5  CONTINUE 
RETURN 
END 
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SUBROUTINE  LEV( C , GAM ,N , AA , LAM , SIGMA) 

SUBROUTINE  TO  CALCULATE  THE  PREDICTOR  FILTER  COEFFICIENTS 
USING  THE  LEVINSON  RECURSION. 

REAL  AA( 0:  100,0:  100) ,C(0: N+2) ,GAM(N) ,LAM( 100) 

REAL  LAMK,SIGMA(0: 100) 

INITIAL  CONDITIONS 
AA( 0,0)  =  1. 

SIGMA(O)  =  C(0) 

COMPUTE  LAM(K) ,  RHO(K);  UPDATE  SIGMA(K)  FOR  THE  NEXT 
ITERATION. 

DO  10  K=1,N 
LAMK  =  0. 

AA(  K ,  0 )  =  1.  0 
DO  20  1=0 ,K-1 

LAMK  =  LAMK  -  C(K-I)*AA(K-1 , I) 

LAM(K)  =  LAMK 
CONTINUE 

GAM(K)  =  LAM( K ) / S I GMA( K - 1 ) 

SIGMA(K)  =  SIGMA(K-l)  -  LAM( K)’VGAM( K) 

COMPUTE  THE  PREDICTOR  FILTER  COEFFICIENTS 
IF(K  .EQ.  1)THEN 
AA( 1,1)  =  GAM(K) 

ELSE 

DO  30  1=1, K-l 

AA( K , I )  =  AA( K- 1 , I )  +  GAM(K)*AA(K-1,K-I) 

CONTINUE 
AA(K,K)  =  GAM(K) 

ENDIF 

CONTINUE 

RETURN- 

END 
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APPENDIX  C.  SPLIT  LATTICE  ALGORITHMS 

PROGRAM  TO  CALCULATE  THE  NTH  ORDER  LATTICE  REFLECTION 
COEFFICIENTS  FROM  A  GIVEN  SEQUENCE  USING  THE  SYMMETRIC  ERROR 
VECTOR,  THE  ASYMMTRIC  ERROR  VECTOR,  OR  THE  FORWARD  AND 
BACKWARD  ERROR  VECTORS.  VARIABLES  DEFINED  IN  PREVIOUS  APPENDICES 
ARE  NOT  REDEFINED. 

*****  THIS  PROGRAM  REQUIRES  SUBROUTINE  INPUT  FROM  APPENDIX  A  ******* 

VARIABLE  DEFINITIONS 

SIG  -  N-TH  DEGREE  NORM  OF  THE  FILTER. 

GAM  -  VECTOR  OF  REFLECTION  COEFFICIENTS  CALCULATED  BY 
THE  LEVINSON  RECURSION. 

LAM  -  REAL  VARIABLE  USED  WHEN  CALCULATING  THE  REFLECTION 
COEFFICIENTS  FROM  THE  LEVINSON  RECURSION 
THE  REFLECTION  COEFFICIENT  IN  TERMS 
OF  THE  FILTER  NORM  IS  GIVEN  BY: 

RHO(K)  =  LAM/SIG 

TAU  -  REAL  VECTOR  OF  "MODIFIED  NORM  VALUES".  THE 
VALUES  ARE  CALCULATED  FROM  A  SUMMATION  OF 
PRODUCT  TERMS  OF  THE  SYMMETRIC  OR  ASYMMETRIC 
PREDICTION  ERROR  SEQUENCES. 

A  -  ARRAY  OF  PREDICTOR  POLYNOMIAL  COEFFICIENTS. 

AS  EACH  I-TH  ROW  OF  THE  ARRAY  CONTAINS  THE 

AL  PREDICTOR  POLYNOMIAL  COEFFICIENTS  FOR  THE  I-TH 

ORDER  PREDICTOR  POLYNOMIAL. 

AR  -  VECTOR  OF  COEFFICIENTS  FROM  THE  KNOWN  TEST  FILTER. 

XO  -  SYMMETRIC  OR  ASSYMMETRIC  PREDICTION  ERROR  VECTOR 
FOR  THE  ( K- 1 )  STAGE  OF  THE  LATTICE  FILTER. 

LL  -  DESIRED  LATTICE  FILTER  ORDER. 

XI  -  SYMMETRIC  OR  ASYMMETRIC  PREDICTION  ERROR  VECTOR 
FOR  THE  K-TH  STAGE  OF  THE  LATTICE  FILTER. 

AT  -  TEMP  STORAGE  FOR  THE  PREDICTION  ERROR  VECTOR  WHILE 
COMPUTING  THE  (K+l)  STAGE  PREDICTION  ERROR  VECTOR. 

FT  -  SHIFTED  FORWARD  PREDICTION  ERROR  VECTOR. 

BT  -  SHIFTED  BACKWARD  PREDICTION  EROR  VECTOR. 

M  -  DESIRED  ORDER  OF  THE  PREDICTOR  POLYNOMIAL. 

T  -  INTEGER  VARIABLE  USED  IN  THE  PROGRAM. 

W  -  WHITE  NOISE  SEQUENCE  VECTOR. 

S  -  INPUT  SEQUENCE  VECTOR 

F  -  FORWARD  PREDICTION  ERROR  VECTOR. 

B  -  BACKWARD  PREDICTION  ERROR  VECTOR. 

VARIABLE  DECLARATIONS 
REAL  AR( 30) ,W(0: 5000), S(0: 5000), RHO( 100) 

REAL  A( 0:  100,0:  100) ,GAM(20) ,RH0S( 100) ,AS(0: 100,0:  100) 

REAL  ALPHA, X1(0: 5000) ,XO(0: 5000) ,AT(0: 5000) ,AL(0: 100,0: 100) 
INTEGER  M,LL, IX,T,L,N 
OPEN( UNIT=4 , BLANK= ' ZERO ' ) 

INITIALIZE  FILTER  ORDER 
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READ(4, 100)M 

C  INITIAL  CONDITIONS  FOR  INPUT  SEQUENCE  GENERATOR 
LL  =  M 
N  =  5000 
IX  =  1913 
WRITE( 6 , 200) 

WRITE( 6 , 300)M 
WRITE! 6 ,400)N 
WRITE! 6 ,500) 

WRITE! 6, 600) 

C  LOOP  TO  GENERATE  WHITE  NOISE  SEQUENCE  AND  TO  READ  TEST  COEFFICIENTS. 
DO  1  1=0 ,N 

CALL  GAUSS! IX, 1.  ,0.  ,V) 

W(I)  =  V 
1  CONTINUE 

C  CALL  STATEMENT  TO  GENERATE  INPUT  SEQUENCE 
CALL  INPUT! LL, W, AR , S ,N) 

C  CALL  STATEMENTS  FOR  SYMMETRIC,  ASYMMETRIC,  AND  LEVINSON  LATTICE 

C  RECURSION  SUBROUTINES 

CALL  SLAT! S,M,N,RHO, ALPHA, XI, AT, XO) 

CALL  AS LAT( S , M , N , RHOS , ALPHA , X 1 , AT , XO ) 

CALL  LEV(N,GAM,M,S) 

WRITE! 6, 700) 

WRITE! 6 , 800 ) 

DO  80  K=1 ,M 

WRITE! 6 ,900)K,RHO(K) ,RHOS!K) ,GAM! K) 

80  CONTINUE 

DO  90  K=1,M 
A(K,0)  =  1. 

AS(K,0)=1. 

AL(K,0)=1. 

IF! K  -EQ.  1)THEN 
A!  1,1)  =  RHO!K) 

AS!  1,1)  =  RHOS( K) 

AL( 1,1)  =  GAM(K) 

ELSE 

DO  95  1=1, K-l 

A(K , I )  =  A(K-1,I)  +  RHO(K)*A(K-l,K-I) 

AS! K  > I )  =  AS!K-1,I)  +  RHOS(K)*AS(K-l,K-I) 

AL(K, I)  =  AL(K-l.I)  +  GAM(K)*AL(K-1,K-I) 

95  CONTINUE 

A(K,K)  =  RHO!K) 

AS!K,K)  =  RHOS(K) 

AL(K,K)  =  GAM!K) 

ENDIF 

90  CONTINUE 

WRITE! 6, 1000) 

WRITE! 6 , 1100) 

WRITE! 6, 1200) 

DO  96  K=1,M 
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WRITE(6 , 1300)M,K,AL(M,K) ,A(M,K) ,AS(M,K) 

96  CONTINUE 

100  FORMAT( 12) 

200  FORMAT(’l') 

300  FORMATC'  FILTER  ORDER  =  ',13) 

400  FORMATC'  V  NUMBER  OF  SAMPLE  POINTS  =  ’,15) 

500  FORMATC ,1 OX, 'KNOWN  FILTER  COEFFICIENTS') 

600  FORMATC '-' ,8X, ' INDEX' ,10X, 'FILTER  COEFFICIENTS') 

700  FORMATC '-' ,10X,' REFLECTION  COEFFICIENTS’) 

800  FORMATC 5X , '  INDEX  ' ,3X, ' SYMMETRIC' ,9X, 'ANTISYMMETRIC' 
+,9X, '  LEVINSON  ') 

900  FORMATC’-’ ,6X, 13 ,6X,F8. 4, 12X.F8.  4, 12X.F8.  4) 

1000  FORMATC '-' ,15X, ’  FILTER  COEFFICIENTS  ') 

1100  FORMATC'-' ,20X,'  LEVINSON  ',12X,'  SPLIT-LEVINSON  ') 

1200  FORMATC 5X, '  INDEX  ’,26X,'  SYMMETRIC  ’ ,4X, '  ASYMMETRIC  ') 

1300  FORMATC'  ' ,6X,2I2, 10X.F8.  4, 11X.F8.  4,7X,F8.  4) 

WRITE( 6 , 200) 

END 

SUBROUTINE  S LAT ( S , M , N , RHO , ALPHA , X 1 , AT , XO ) 

C  SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 

C  USING  THE  SYMMETRIC  PREDICTION  ERROR  VECTOR. 

REAL  X 1 ( 0 : M+N) ,XOCO:  M+N) ,RHO(M) ,S(0:  N) .ALPHA 
REAL  AT(0:  M+N), A( 100,100) 

C  INITIAL  CONDITIONS 
INTEGER  T 
RRHO  =  0. 

TAU  *  0. 

C  INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C  STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0 ,N- 1 
XOCT)  =  2.*S(T) 

TAU  =  TAU  +  SCT)**2 
10  CONTINUE 

DO  20  T=0,N 
IFCT. EQ. N)SCT)  =  0. 

IF(T.  EQ.  0)THEN 
X1(T)  =  SCT) 

ELSE 

X1CT)  =  SCT)  +  S(T-l) 

ENDIF 

20  CONTINUE 

C  LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  101  K=1,M 
TTAU  =  TAU 

C  STORE  TAU(K-l) ,  AND  COMPUTE  TAU(K) . 

TAU  =  0. 

DO  30  T=0 ,N+K-1 
TAU  =  TAU  +  X1(T)**2 
30  CONTINUE 
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TAU  =  TAU/2. 

C  COMPUTE  ALPHA(K) ,  RHO(K);  STORE  TAU(K)  AND  RHO(K)  FOR 
C  NEXT  ITERATION. 

ALPHA  =  TAU/TTAU 
TTAU  =  TAU 

RHO(K)  =  1.  -  (ALPHA/( 1.  +  RRHO) ) 

RRHO  =  RHO(K) 

C  LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 

DO  40  T=0 ,  N+K 
IF(T  . EQ.  0)THEN 
AT(T)  =  X1(T) 

ELSEIF(T. EQ.  N+K)THEN 
AT(T)  =  Xl(T-l) 

ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA*XO(T-l) 

END  IF 

40  CONTINUE 

C  LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 
C  (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 

DO  50  T=0 , N+K- 1 
XO(T)  =  X1(T) 

50  CONTINUE 

DO  60  T=0 ,N+K 
Xl( T)  =  AT(T) 

60  CONTINUE 

101  CONTINUE 

RETURN 
END 

SUBROUTINE  ASLAT( S , M , N , RHOS , ALPHA , XI , AT , XO ) 

C  SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 

C  USING  THE  ASYMMETRIC  PREDICTION  ERROR  VECTOR. 

REAL  Xl( 0: M+N) ,XO(0: M+N) ,RHOS(M) ,S( 0: N) .ALPHA 
REAL  AT( 0 : M+N ) 

INTEGER  T 

C  INITIAL  CONDITIONS 
RRHO  =  0. 

TAU  =  0. 

C  INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C  STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0,N-I 
XO(T)  =  0. 

TAU  =  TAU  +  S(T)**2 
10  CONTINUE 

DO  20  T=0,N 

IF(T. EQ. N)X1(T)  =  -S(T-l) 

IF(T. EQ. 0)THEN 
X1(T)  =  S(T) 

ELSE 

X1(T)  =  S(T)  -  S(T-l) 

ENDIF 
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20  CONTINUE 

C  LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  101  K=1,M 

C  STORE  TAU(K-l),  AND  COMPUTE  TAU(K). 

TTAU  =  TAU 
TAU  =  0. 

DO  30  T=0 ,N+K- 1 
TAU  =  TAU  4-  X1(T)**2 
30  CONTINUE 

TAU  =  TAU/ 2. 

C  COMPUTE  ALPHA(K) ,  RHO(K);  STORE  TAU(K)  AND  RHO(K)  FOR 
C  NEXT  ITERATION. 

ALPHA  =  TAU/TTAU 
TTAU  =  TAU 

RHOS(K)  =  (ALPHA/C  1.  -  RRHO) )  -  1. 

RRHO  =  RHOS(K) 

C  LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 

DO  40  T=0,N+K 
IF(T  . EQ.  0)THEN 
AT(T)  =  X1(T) 

ELSEIFCT.  EQ.  N+K)THEN 
AT(T)  =  Xl(T-l) 

ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA*XO(T-l) 

END  IF 

40  CONTINUE 

C  LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 

C  (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 

DO  50  T=0,N+K-1 
XO(T)  *  X1(T) 

50  CONTINUE 

DO  60  T=0 , N+K 
X1(T)  =  AT(T) 

60  CONTINUE 

101  CONTINUE 

RETURN 
END 

SUBROUTINE  LEV(N,GAM,M,S) 

C  SUBROUTINE  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS  FOR  AN 

N-TH  ORDER  LATTICE  FILTER  FROM  THE  FORWARD  AND  BACKWARD 

PREDICTION  ERROR  VECTORS. 

REAL  F( 0:  5100) ,B(0:  5100) ,FT(0:  5100) , BTC  0:  5100) ,GAM(20) 

REAL  LAM,SIG,S(0:  N) 

INTEGER  T 

C  INITIAL  CONDITIONS;  INITIALIZE  FORWARD  AND  BACKWARD  PREDICTION 

C  ERROR  VECTORS. 

SIG  =  0. 

DO  10  1=0, N-l 
F(I)  =S(I) 

B(I)  =  S(I) 
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FT(I)  =  S(I) 

SIG  =  SIG  +  S(I)**2 
10  CONTINUE 

C  LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  20  K=1,M 

C  FORM  SHIFTED  ERROR  VECTORS 
DO  30  T=1 , N+K- 1 
BT(T)  =  B(T-l) 

30  CONTINUE 

BT( 0)  =  0. 

FT( N+K- 1 )  =  0. 

LAM  =  0. 

COMPUTE  LAM(K) ,  GAM(K);  UPDATE  K-TH  ERROR  NORM  AND 
STORE  FOR  NEXT  ITERATION. 

DO  40  T=l,N+K-2 
LAM  =  LAM  -  FT(T)*BT(T) 

40  CONTINUE 

GAM(K)  =  LAM/SIG 

IF( K  .  EQ.  M)GOTO  20 
SIG  =  SIG  -  LAM*GAM(K) 

C  COMPUTE  ( K+l)  FORWARD  AND  BACKWARD  PREDICTION  ERRORS  AND  SHIFT 
C  INTO  TEMPORARY  VECTORS  FOR  NEXT  ITERATION. 

DO  50  T=0 ,N+K- 1 
F(T)  =  FT(T)  +  GAM(K)*BT(T) 

B(T)  =  GAM( K)*FT(T)  +  BT(T) 

50  CONTINUE 

DO  60  T=0 ,N+K-1 
FT(T)  =  F(T) 

BT(T)  =  B(T) 

60  CONTINUE 

20  CONTINUE 

RETURN- 
END 
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APPENDIX  D.  MA  PREDICTOR  COEFFICIENT  PROGRAM 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  IS  TO  COMPUTE  THE  FIR  FILTER  COEFFICIENTS 
USING  THE  SPLIT-LEVINSON  ALGORITHM,  AND  AN  AUTOREGRESSIVE 
MOVING  AVERAGE  PROCESS.  VARIABLES  DEFINED  IN  PREVIOUS 
APPENDICES  ARE  NOT  REDEFINED. 

****  this  PROGRAM  REQUIRES  THE  SUBROUTINES  ODD,  AND  EVEN 
FROM  APPENDIX  A  **** 

VARIABLE  DEFINITIONS 

ENORM  -  PREDICTOR  COEFFICIENT  ERROR  NORM. 

ENORM  =  SQRT( ( A( 0) -AA( 0) )**2  +. . . +  (A(N)-AA(N) )**2) 
NSTRT  -  NUMBER  OF  POINTS  OF  INPUT  SEQUENCE  TO  START. 

NEND  -  NUMBER  OF  POINTS  OF  INPUT  SEQUENCE  AT  END  OF  PROGRAM. 

DELX  -  ERROR  VECTOR. 

DELY  -  ERROR  VECTOR. 

NPTS  -  NUMBER  OF  INPUT  DATA  POINTS  (0,1,... ,NPTS). 

EXB  -  BACKWARD  PREDICTION  ERROR  POWER  OF  X. 

RXY  -  VECTOR  OF  X  AND  Y  CROSSCORRLATION  ELEMENTS 

AA  -  VECTOR  OF  CALCULATED  PREDICTOR  COEFFICIENTS. 

NX  -  INDEX  FOR  X-AXIS  OF  PLOTTING  ROUTINE. 

EN  -  VECTOR  OF  PREDICTOR  COEFFICIENT  NORMS. 

EX  -  FORWARD  PREDICTIZN  ERROR  POWER  OF  X. 

EY  -  FORWARD  PREDICTION  ERROR  POWER  OF  Y. 

KY  -  Y  REFLECTION  COEFFICIENT. 

KX  -  X  REFLECTION  COEFFICIENT. 

LL  -  FILTER  ORDER  VARIABLE  USED  IN  SUBROUTINE  CORR. 

IX  -  INTEGER  SEED  NUMBER  USED  BY  IMSL  SUBROUTINE  GAUSS. 

RX  -  X  AUTOCORRELATION  VECTOR. 

RY  -  Y  AUTOCORRELATION  VECTOR. 

B  -  MA  COEFFICIENT  VECTOR. 

C  -  MA  COEFFICIENT  VECTOR. 

D  -  MA  COEFFICIENT  VECTOR. 

T  -  INTEGER  VARIABLE  USED  IN  THE  SUBROUTINE  ODD 

IN  THE  COMPUTATION  OF  THE  TAU(K)'S. 

L  -  DUMMY  VARIABLE  USED  DURING  THE  CALCULATION 

OF  THE  SYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIAL 
COEFFICIENTS. 

X  -  INPUT  WHITE  NOISE  VECTOR. 

M  -  INDEXING  VARIABLE  USED  IN  FIR  FILTER  COEFFICIENT 

RECURSION  (M=l ,2 , . . . ,N). 

Y  -  OUTPUT  SEQUENCE  FROM  FIR  TEST  FILTER. 

VARIABLE  DECLARATIONS 

REAL  P( 0:  100,0:  100) ,TAU(0:  100) ,C(0: 50) ,3(0: 50) ,EN(200) 

REAL  A( 0:  100,0: 100) ,LAMDA(0:  100) ,X(0:  10000) ,D(0:  50) 

REAL  AR( 0: 30) ,Y(0: 10000) ,EY(0: 50) ,EX(0: 50) ,KX(50) 

REAL  DELX( 0:  50) ,DELY( 0:  50) ,EXB(0: 50) ,KY(50) ,XX(200) 

REAL  RX( 0: 100) , AL?HA( 100) ,RY(0: 2) ,RXY(0: 100) ,AA(0: 100) 
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INTEGER  M , LL , IX ,T , KODD , KEVES , L , N , NPTS 
C  DESIRED  FILTER  ORDER  AND  THE  TEST  FILTER  COEFFICIENTS 
C  ARE  READ  FROM  A  DATA  FILE  (FILE  FT04F001). 

0PEN(lTNIT=4 , BLANK=' ZERO’  I 

C  INITIALIZE  FILTER  ORDER  AND  TEST  FILTER  COEFFICIENTS 
READ(4, 100)N 
LL  =  N 

DO  20  K=0 , LL 
READ( 4 , 700 ) AR( K) 

20  CONTINUE 

NPTS  =  4000 
WRITE( 6 , 300) 

WRITE! 6 , 400 )N 
IX  =  1913 
WRITE! 6 , 600 ) NPTS 

C  GENERATE  (NPTS+1)  WHITE  NOISE  SEQUENCE 
DO  10  1=0, NPTS 
CALL  GAUSS( IX, 1.  ,0.  ,V) 

X(I)  =  V 
10  CONTINUE 

C  CREATE  TNP'T  SEQUENCE  FROM  GIVEN  FIR  TEST  FILTER 
CALL  INPUT! LL , X , AR , Y ,NPTS) 

C  INITIALIZE  AUTOCORRELATION  VECTOR  RX.RY,  AND  CROSSCORRELATION 
C  VECTOR  RXY 

CALL  CORR( N+l , NPTS , X , Y , RX , RY , RXY ) 

WRITE( 6 , 800) 

DO  30  1=0, N 

WRITE! 6,900)1, AR(I) ,RX( I ) 

30  CONTINUE 

C  INITIAL  CONDITIONS  FOR  MOVING  AVERAGE  MODEL  PARAMETERS 
BOO  =  -RXY(0)/RX(0) 

EY( 0 )  =  RY(0)  -  B00*RXY( 0) 

EX( 0 )  =  RX( 0) 

DO  40  1=0,1 
C(I)  =  I 
D(  I)  =  I 
40  CONTINUE 

B(0)  =  1. 

B ( 1 )  =  BOO 

DELYCO)  =  RXY(l)  -  B00*RX(1) 

DELX(O)  =  RX( 1) 

C  LOOP  TO  GENERATE  PREDICTOR  COEFFICIENT  VECTOR 
DO  120  M=1,N 

C  INITIAL  CONDITIONS  FOR  THE  SYMMETRIC  PREDICTOR  POLYNOMIAL 
C  CALCULATIONS. 

P(0,0)  =  2. 

P(l,l)  =  I- 

TAU(O)  =  RX(0) 

LAMDA(O)  =  1. 

KODD  =  1 
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KEVEN  =  2 
A(M,0)  =  1. 

C  SYMMETRIC  SPLIT-LEVINSON  RECURSION 
DO  70  K=1,M 
P(K,0)  =  1. 

C  CALL  STATEMENTS  FOR  EVEN  OR  ODD  VALUES  OF  INDEX  K 
IF( K. EQ.  KODD)THEN 

CALL  ODD( RX , P , K , N , TAU , KODD , T ) 

ELSEIF(K. EQ. KEVEN) THEN 
CALL  EVENT  RX ,  P ,  K ,  N , TAU , KEVEN , T) 

END  IF 

ALPHA(K)  =  TAUT  K)/TAU( K- 1 ) 

C  LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  60  1=1, T 

P( K+l , I )  =  P( K , I )  +  PCK.I-l)  -  ALPHA(K)*P( K- I ,1-1) 

C  DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IF(I.LT.  T  .OR.  I.EQ.  K)GOTO  60 
L  =  K+l 
DO  50  J=I+1 ,K 
P( L , J )  =  P(L,L-J) 

50  CONTINUE 

60  CONTINUE 

LAMDA(K)  =  2.  -  ( ALPHA! K)/LAMDA(K-1) ) 

70  CONTINUE 

C  CALCULATION  OF  N-TH  ORDER  PREDICTOR  COEFFICIENTS ,  A(K, 1 ) , .  .  .  A(K,K) 

C  COMPUTE  PREDICTOR  COEFFICIENTS  FOR  K-TH  ITERATION 

DO  80  1=1, M 

A(M, I)  =  A(M, 1-1)  +  P(M+1 , I)  -  LAMDA(M)*P(M,I-1) 

C( I )  =  ACM, I) 

80  CONTINUE 

DO  90  J=1,M 
D( 1)  =  -C(J) 

IF(  J  .LT.  M)THEN 
DO  95  I=J+1,M 
DC  I)  »  DC  I -1) 

95  CONTINUE 

ENDIF 

90  CONTINUE 

D( 0)  =  0. 

D(  M+l )  =  1. 

C  UPDATE  PREDICTION  ERRORS 
XBTMP  =  0. 

XTMP  =  0. 

DO  25  1=1, M 

XBTMP  =  XBTMP  +  C( I)*RXY(M+1-I) 

XTMP  =  XTMP  +  C(I)*RX(M+1-I) 

25  CONTINUE 

DELX(M)  =  RX(M+1)  -  XTMP 
EXB(M)  =  RXY( M+l )  -  XBTMP 
C  UPDATE  REFLECTION  COEFFICIENTS 
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KX( M)  =  -DELX(M-1)/EX(M-1) 

EX(  M)  =  EX(M-l)  +  KX ( M ) *DELX( M - 1 ) 

KY(M)  =  -DELY( M- 1 ) / EX( M) 

EY(M)  =  EY(M-l)  +  KY(>1)*EXB(M) 

C  UPDATE  B  VECTOR 
B( M+l )  =  0. 

DO  45  1=0, M+l 
B(I)  =  B(.I)  +  KY(M)*D(D 
45  CONTINUE 

YTMP  =  0. 

DO  55  1=1, M+l 

YTMP  =  YTMP  -  B( I)*RX(M+2-I) 

55  CONTINUE 

DELY(M)  =  RXYCM+l)  -  YTMP 
IF(M  .EQ.  N)THEN 
WRITE( 6 , 1000 )N 
WRITE( 6 , 1100 ) 

DO  130  K=0 , M 
WRITE( 6 , 1200)K, -B(K+1) 

AA(K)  =  - B ( K+ 1 ) 

130  CONTINUE 
END  IF 

120  CONTINUE 

100  FORMATC 12) 

200  FORMATC '  ' ) 

300  FORMAT('l') 

400  F0RMAT('  FILTER  ORDER  =' ,13) 

500  FORMATC'-') 

600  FORMATC’-’,’  NUMBER  OF  SAMPLE  POINTS  =  ’,15) 

700  FORMATC F8. 4) 

SCO  FORMATC '5X, ’ INDEX' ,5X, 'KNOWN  COEFFICIENTS ', 5X , 

+’ AUTOCORRELATION  FUNCTION  RCK) ' ) 

900  FORMATC'  ' , 5X , 13 , 1 IX ,F8. 4 , 2 IX ,F8. 4) 

1000  FORMATC ,10X, 'PREDICTOR  COEFFICIENTS  FOR  FILTER  ORDER  =  ’,13) 
1100  FORMATC ,5X, ' INDEX’ , 12X, 'FIR  PREDICTOR  COEFFICIENTS') 

1200  FORMATC'  ' ,5X, 13 ,23X,F8. 4) 

WRITEC  6 , 300) 

END 

SUBROUTINE  INPUT(LL,X,AR,Y,NPTS) 

SUBROUTINE  TO  GENERATE  THE  INPUT  SEQUENCE  FROM  A  GIVEN  FIR 
FILTER  AND  ZERO  MEAN,  UNIT  VARIANCE  WHITE  NOISE. 

REAL  XC 0:  NPTS) ,AR(0:  LL) ,Y(0: NPTS) 

C  CALCULATE  INPUT  SEQUENCE  VALUES  BETWEEN  0  AND  FILTER  ORDER. 

DO  45  K=0 , NPTS 
TEMP  =  0. 

IFCK. LE. LL)THEN 
LU=K 
ELSE 
LU=LL 
ENDIF 
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J  =  K 

DO  55  1=0, LU 

TEMP  =  TEMP  +  AR(I)*X(J) 

J  =  J-l 
55  CONTINUE 
Y(K)  =  TEMP 
45  CONTINUE 
RETURN 
END 

SUBROUTINE  CORR( NLAG , NPTS , X , Y , RX , RY , RXY) 

C  SUBROUTINE  TO  CALCULATE  THE  AUTOCORRELATION  FUNCTION  X  AND  Y  AND 

C  THE  CROSSCORRELATION  FUNCTION  BETWEEN  X  AND  Y 

REAL  RX(  0:  NLAG) ,Y(0:  NPTS) ,X(0:  NPTS) ,RXY(0:  NLAG) ,RY(0:  2) 
INTEGER  T 

C  COMPUTE  THE  AUTOCORRELATION  OF  X  AND  THE  CROSSCORRELATION  OF 
C  X  AND  Y  FOR  LAGS  0,1,2,... ,NLAG 
DO  5  1=0, NLAG 
XTEMP  =  0. 

XYTEMP  =  0. 

C  COMPUTE  THE  AUTOCORRELATION  OF  X  AND  THE  CROSSCORRELATION  OF 
C  X  &  Y  FOR  LAG  I 

DO  15  T=0,NPTS-1-I 
XTEMP  =  XTEMP  +  X(T)*X(T+I) 

XYTEMP  =  XYTEMP  +  X(T)*Y(T+I) 

15  CONTINUE 

RX( I )  =  XTEMP* ( 1.  /FLOAT(NPTS-l-I) ) 

RXY ( I )  =  XYTEMP* ( 1. /FLOAT(NPTS-l-I) ) 

C  COMPUTE  THE  ZERO  LAG  AUTOCORRELATION  FUNCTION  OF  Y 
IF( I  . EQ.  0)THEN 
RY( 0)  =  0. 

DO  16  J=0 , NPTS- 1 
RY( 0 )  =  RY( 0)  +  Y(J)**2 

16  CONTINUE 

RY( 0 )  =  RY(0)*( 1. /FLOAT(NPTS-l) ) 

END  IF 

5  CONTINUE 

RETURN 
END 
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APPENDIX  E.  EXTENDED  PROXY  PROGRAM 

PROGRAM  TO  CALCULATE  THE  NTH  ORDER  LATTICE  REFLECTION 
COEFFICIENTS  FROM  A  GIVEN  SEQUENCE  USING  THE  SYMMETRIC  ERROR 
VECTOR,  THE  ASYMMTRIC  ERROR  VECTOR,  OR  THE  FORWARD  AND 
BACKWARD  ERROR  VECTORS. 

VARIABLE  DEFINITION'S 

PT  -  TEMPORARY  ARRAY  USED  TO  AVERAGE  PREDICTOR 
COEFFICIENTS. 

PP  -  ESTIMATED  NUMBER  OF  COMPLEX  SINUSOIDS  PRESENT. 

A1  -  AMPLITUDE  OF  ifl  SINUSOID,  (1-4)  SINUSOIDS 
PRESENT. 

FS  -  SAMPLING  FREQUENCY. 

FI  -  FREQUENCY  OF  #1  SINUSOID  IN  TEST  SEQUENCE. 

THETA 1  -  DIGITAL  FREQUENCY  OF  #1  TEST  ANALOG  FREQUENCY. 

VARIABLE  DECLARATIONS 

REAL  W( O:  5000) ,S(0:  5000) , ALPHA( 100) ,ROOTR( 100) ,XCOF(0:  100) 

REAL  PC  0:  100,0:  100) ,ALPHAS( 100) ,COF( 0:  100) ,ROOTI( 100) 

REAL  X 1 ( 0:  5000), X0(0:  5000) ,AT(0: 5000) ,PS(0:  100,0:  100) ,PT(0:  100) 
INTEGER  T , PP 

OPENC UNIT=4 , B  LANK= ’ZERO') 

INITIALIZE  FILTER  ORDER 
READ(4, 1Q0)PP 
M  =  2*PP 

INITIAL  CONDITIONS  FOR  INPUT  SEQUENCE  GENERATOR 
IX  =  1913 
A  =  SQRTC2. ) 

B  =  SQRT( 2. ) 

C  =  SQRT( 10.  ) 

D  =  SQRT( 2.  0) 

E  =  SQRT(2.0) 

FI  =  5. 5E+01 
F2  =  7. 5E+02 
F3  =  1. 25E+02 
F4  =  1. 75E+02 
FS  =  2. 25E+02 
TWOPI  =  6. 2831853 
THETA1  =  (TW0PI*F1 )/FS 
THETA2  =  (TVOPI*F2)/FS 

LOOP  TO  GENERATE  WHITE  NOISE  SEQUENCE  AND  TO  READ  TEST  COEFFICIENT 
DO  1  1=0, N 

CALL  GAUSS( IX , 1. ,0. ,V) 

W( I )  =  c*v 

A1  =  A*COS(TWOPI*(Fl/FS)*FLOAT(I)) 

A2  =  B*COS(TWOPI*(F2/FS)*FLOAT(I)) 

A3  =  D*COS(TWOPI*(F3/FS)*FLOAT(I)) 

A4  =  E*C0S(TW0PI*(F4/FS)*FL0AT(I)) 

S(I)  =  A1  +  A2  +  W(I) 
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I  CONTINUE 

C  CALL  STATEMENTS  FOR  SYMMETRIC,  ASYMMETRIC,  AND  LEVINSON  LATTICE 
C  RECURSION  SUBROUTINES 

CALL  SLAT( S ,M , N , P , ALPHA , XI , AT , XO ) 

C  CALL  ASLAT( S ,M,N, PS .ALPHAS , XI , AT,XO) 

WRITEC 6 , 200) 

WRITEC  6 , 300 )PP 
WRITE(6,400)M 
WRITE( 6 , 500)N 

C  DISPLAY  COEFFICIENTS  OF  POLYNOMIAL 

.RITE(6,600) 

DO  11  K=0 , M 

IF( K  .  EQ.  M)P(M,K)  =  1. 0 

WRITEC  6,700)M,K,P(M,K) 

II  CONTINUE 

100  FORMATC 14) 

200  FORMATC  1') 

300  FORMATC'  NUMBER  OF  COMPLEX  EXPONENTIALS  IN  MODEL  =  ’,13) 
400  FORMATC'  SYMMETRIC  FILTER  ORDER  =  ’,13) 

500  FORMATC'  V  NUMBER  OF  SAMPLE  POINTS  =  ’,15) 

600  FORMATC'-' ,8X, 'INDEX' , 13X, ’ COEFFICIENTS' ) 

700  FORMATC 5X, 212, 16X,F8. 4) 

WRITEC 6, 200) 

END 

SUBROUTINE  SLATCS ,M,N,P, ALPHA, XI ,AT,XO) 

C  SUBROUTINE  TO  COMPUTE  THE  SYMMETRIC  PREDICTOR  COEFFICIENTS 

C  USING  THE  SYMMETRIC  PREDICTION  ERROR  VECTOR. 

REAL  X 1 ( 0:  M+N) ,XOCO:  M+N) , ALPHAC M) , S( 0: N) 

REAL  AT( 0: M+N), P(0: 100,0: 100) 

INTEGER  T 

C  INITIAL  CONDITIONS 

KODD  =  1 
KEVEN  =  2 
TAU  =  0. 

P(  1 , 1)  =  1. 

P(0,0)  =  2. 

C  INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C  STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0 , N -  1 
XO(T)  =  2.*S(T) 

TAU  =  TAU  +  S ( T ) **  2 
10  CONTINUE 

DO  20  T=0,N 
IF(T. EQ. N)S(T)  =  0. 

IF(T. EQ.  0)THEN 
X1(T)  =  S(T) 

ELSE 

X1(T)  =  S(T)  +  SCT-1) 

END  IF 
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CONTINUE 

C  LOOP  TO  COMPUTE  THE  SYMMETRIC  PREDICTOR  COEFFICIENTS 
DO  101  K=1 ,M 
P(K,0)  =  1.  0 
TTAU  =  TAU 
IF(K  . EQ.  KODD)THEN 
LL  =  (K+l)/2 
ELSE 
LL  =  K/2 
END  IF 

C  STORE  TAU(K-l)  ,  AND  COMPUTE  TAU'(K). 

TAU  =  0. 

DO  30  T=0 ,N+K- 1 
TAU  =  TAU  +  X1(T)**2 
30  CONTINUE 

TAU  =  TAU/ 2. 

C  COMPUTE  ALPHA(K);  STORE  TAU(K)  FOR  NEXT  ITERATION. 

ALPHA(K)  =  TAU/TTAU 
TTAU  =  TAU 

C  LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  40  1=1, LL 

PCK+1,1)  =  P(K, I)  +  P( K , I  - 1 )  -  ALPHA! K)*P(K- 1 ,1-1) 

C  DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IF(I.  LT.  LL  .  OR.  I.  EQ.  K)GOTO  40 
L  =  K+l 
DO  50  J=I+1,K 
P(L,J)  =  P(L,L-J) 

50  CONTINUE 

40  CONTINUE 

C  LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 

DO  60  T=0 , N+K 
IF(T  .EQ.  0)THEN 
AT(T)  =  X 1 ( T ) 

ELSEIFCT.  EQ. N+K) THEN 
AT(T)  =  Xl(T-l) 

ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA! K)*XO(T-l) 

END  IF 

60  CONTINUE 

C  LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 

C  (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 

DO  70  T=0 , N+K - 1 
XO(T)  =  X1(T) 

70  CONTINUE 

DO  80  T=0 ,N+K 
X1(T)  =  AT(T) 

80  CONTINUE 

IF(K  .EQ.  KODD)THEN 
KODD  =  KODD  +  2 
ELSE 
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KEVEN  =  KEVEN  +  2 
END  IF 

101  CONTINUE 
RETURN 
END 

SUBROUTINE  ASLAT( S , M , N , PS , ALPHAS , X 1 , AT , XO ) 

C  SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 

C  USING  THE  ASYMMETRIC  PREDICTION  ERROR  VECTOR. 

REAL  X 1 ( 0 : M+N) ,XO(0: M+N ) , PS( 0:  100,0:  100) , SCO: N) ,ALPHAS(M) 
REAL  AT( 0:  M+N) 

INTEGER  T 

C  INITIAL  CONDITIONS 

PS(0,0)  =  2. 

PS (  1,1)  =  1. 

KODD  =  1 
KEVEN  =  2 
TAU  =  0. 

C  INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C  STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0 ,N- 1 
XO(T)  =  0. 

TAU  =  TAU  +  S(T)**2 
10  CONTINUE 

DO  20  T=0 ,N 

IF(T.  EQ.N)Xl(T)  =  -S(T-l) 

IF(T.  EQ.  0)THEN 
X1(T)  =  S(T) 

ELSE 

X1(T)  =  S(T)  -  S(T-l) 

END  IF 

20  CONTINUE 

C  LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  101  K=1 ,M 
PS(K,0)  =  1. 

TTAU  =  TAU 
IF(K  .EQ.  KODD)THEN 
LL  =  (K-1J/2 
ELSE 

LL  =  K/2 
END  IF 

C  STORE  TAU(K-l) ,  AND  COMPUTE  TAU(K). 

TAU  =  0. 

DO  30  T=0 ,N+K- 1 
TAU  =  TAU  +  X1(T)**2 
30  CONTINUE 

TAU  =  TAU/ 2. 

C  COMPUTE  ALPHA(K)  ;  STORE  TAU(K)  FOR  NEXT  ITERATION. 

ALPHAS(K)  =  TAU/TTAU 
TTAU  =  TAU 

C  LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
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DO  40  1=1, LL 

PS( K+l , I )  =  PS( K , I )  +  PS(K.I-l)  -  ALPHAS! K)*PS( K- 1 , I -1 ) 

C  DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IFCI.  IT.  LL  .  OR.  I.EQ. K,OGTC  40 
L  =  K+l 
DO  50  J=I+1,K 
PS(L.J)  =  -PS( L , L- J) 

50  CONTINUE 

40  CONTINUE 

C  LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 

DO  60  T=0 , N+K 
IF(T  .  EQ.  0)THEN 
AT(T)  =  X 1 ( T ) 

ELSEIF(T. EQ. N+K) THEN 
AT( T)  =  Xl(T-l) 

ELSE 

AT(T)  =  Xl( T)  +  Xl(T-l)  -  ALPHAS! K)*XO(T- 1 ) 

ENDIF 

60  CONTINUE 

C  LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 

C  (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 

DO  70  T=0 , N+K- 1 
XO( T)  =  Xl( T) 

70  CONTINUE 

DO  80  T=0 , N+K 
X1(T)  =  AT( T) 

80  CONTINUE 

IF( K  .EQ.  KODD)THEN 
KODD  =  KODD  +  2 
ELSE 

KEVEN  =  KEVEN  +  2 
ENDIF 

101  CONTINUE 

RETURN- 
END 
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